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E.  A.  HUMPHREYS 


Materials  Sciences  Corporation 
SUMMARY 

This  report  describes  the  development  and  implementation  of  a 
methodology  to  predict  damage  initiation  and  growth  in  composite 
laminates  when  subjected  to  low  velocity,  low  mass,  lateral  impact. 
The  methodology  incorporates  a  transient  dynamic  finite  element 
analysis  with  composite  stress  and  failure  analyses.  The  procedure 
incorporates  damage,  as  it  is  predicted,  into  the  displacement  solu¬ 
tion.  Thus,  the  damage  predicted  during  any  time  step  is  incorporated 
into  the  dynamic  solution  at  future  time  steps.  This  coupling  of 
damage  and  dynamic  response  is  the  heart  of  the  computerized 
procedure . 

Utilizing  a  perfectly  plastic  impact  assxamption,  the  impact 
phenomenon  is  reduced  to  an  initial  velocity  dynamics  problem  with 
the  impacting  mass  lumped  at  appropriate  locations.  The  displace¬ 
ment  time  response  of  the  laminated  plate  is  predicted  using  the 
transient  analysis.  Composite  ply  stresses  and  interlaminar  shear 
stresses  are  computed  based  on  nodal  moments  and  forces  and  laminate 
models.  Failure  analyses  are  performed  and  appropriate  elemental 
properties  degraded  within  the  displacement  solution  matrices. 

The  analysis  procedure  has  been  utilized  to  simulate  the  impact 
of  a  1.59  cm.  steel  sphere  on  an  eight  ply,  [45/0/-45/90] g  Gr/Ep 
laminate.  The  damage  predicted  included  interlaminar  shear  fail¬ 
ures,  laminate  back  face  splitting  and  the  progression  of  damage 
within  the  plane  of  the  plate  and  through  the  thickness  of  the  plate. 
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I .  INTRODUCTION 


With  the  ever  increasing  use  of  laminated  composites  in  struc¬ 
tural  applications,  an  interesting  phenomenon  has  become  apparent. 
This  phenomenon  concerns  the  real  possibility  of  invisible  damage 
within  a  composite  structure  caused  by  low  velocity,  low  mass  im¬ 
pacts  .  This  type  of  loading  environment  is  most  easily  envisioned 
as  the  impact  of  a  dropped  workman's  tool  on  a  structure  or,  per¬ 
haps,  runway  debris  ejected  onto  a  structure. 

The  primary  concern  related  to  this  type  of  loading  environ¬ 
ment  is  the  introduction  of  performance  degrading  damage  within  the 
composite  structure.  This  damage  may  not  be  visually  apparent  dur¬ 
ing  subsequent  inspections  of  the  structure  and,  hence,  the  perform¬ 
ance  degradation  will  also  not  be  immediately  apparent. 

The  subject  of  impact  related  phenomena  has  been  studied  by 
many  investigators  utilizing  many  different  approaches.  Much  of 
this  work  has  been  related  to  ballistic  type  impact  and,  hence,  is 
not  applicable  here.  In  ballistic  impact,  the  velocities  involved 
are  high  enough  to  promote  large  stress  wave  propagation  effects . 

Prior  work  in  the  analysis  of  the  low  speed  impact  problem  has 
established  that  it  is  reasonable  to  neglect  the  stress  wave  propa¬ 
gation  problem  and  to  focus  on  the  transient  structural  dynamic  ap¬ 
proach.  Different  approaches  have  appeared  in  the  literature  to 
combine  contact  effects  with  dynamic  effects.  The  Hertz  contact 
problem  has  been  extended  to  the  problem  of  dynamic  contact  and  also 
to  the  problem  of  contact  of  anisotropic  bodies  (see  ref.  1) .  These 
approaches  treat  impact  with  a  semi-infinite  target.  In  the  present 
case,  one  is  concerned  with  a  target  in  which  the  dynamic  response 
of  the  target  is  important  in  the  sense  of  transient  structural 
motion  rather  than  material  displacement.  This  problem  appears  to 
have  been  addressed  first  by  Timoshenko  (ref.  2),  as  described  by 
Goldsmith  (ref.  3) .  Timoshenko  studied  the  problem  of  the  impact  of 
a  beam  where  the  contact  between  the  bodies  was  governed  by  Hertz ' s 
law  for  contact  deformations.  Karas  (ref.  4)  extended  the  Timoshenko 


approach  to  the  study  of  plate  impact  (see  ref .  5) .  Moon  (refs.  6, 

7)  has  utilized  the  Hertz  impact  theory  in  combination  with  a  Mindlin 
plate  theory  (ref.  8)  to  model  a  similar  approach  for  impact  of  plate 
structures.  This  particular  approach  yields  a  nonlinear  mathematical 
problem  and  extended  numerical  analysis  is  required  to  obtain  solu¬ 
tions.  The  procedure  is  sufficiently  complex  to  motivate  considera¬ 
tion  of  alternate  approaches. 

Two  such  approaches  are  based  on  simplifications  of  the  contact 
force  analysis.  In  one  case,  it  is  considered  that  the  impact  takes 
place  during  a  time  period  which  is  very  short  compared  to  the  period 
of  the  first  natural  frequency.  In  this  case,  it  appears  reasonable 
to  regard  the  impact  as  having  imparted  an  impulse  locally  to  the 
plate  and  to  then  study  the  dynamic  response  of  the  composite  plate 
target  to  that  impulsive  loading.  This  approach  (see  ref.  9)  is 
appropriate  as  the  structural  stiffness  increases  and  the  impacting 
mass  decreases. 

Another  line  of  approach  initiated  by  Clebsch  (ref.  10)  as  de¬ 
scribed  by  Goldsmith  (ref.  3)  assumes  that  upon  impact,  the  projec¬ 
tile  moves  with  the  plate  and  that  the  velocity  of  the  projectile 
becomes  an  initial  velocity  condition.  Thus,  the  analysis  is  the 
structural  dynamic  response  of  the  plate  with  the  attached  mass. 
McQuillen  et  al.  (ref.  11)  applied  this  approach  to  a  beam.  They 
minimized  some  of  the  nximerical  problems  by  considering  the  contact 
zone  between  projectile  and  target  to  have  finite  width.  This  ap¬ 
proach  tended  to  minimize  the  contribution  of  the  higher  frequency 
modes  and,  thus,  numerical  procedures  were  more  successful.  How¬ 
ever,  even  with  these  assumptions,  the  work  of  reference  11  shows 
that  modes  of  vibration  other  than  the  fundamental  mode  can  be  ex¬ 
cited  by  impact,  particularly  if  the  striker  mass  is  small.  This 
approach,  which  is  expanded  somewhat  in  references  12  and  13,  is 
being  utilized  in  the  current  effort. 

In  the  present  study,  the  primary  emphasis  has  focused  on  the 
prediction  of  damage  initiation  and  propagation  during  the  impact 
event  and  subsequent  dynamic  plate  response.  The  difficulty  posed 


here  is  that  any  induced  damage  will  alter  the  plate  stiffness  lo¬ 
cally  and,  hence,  affect  the  subsequent  dynamic  response.  Thus, 
the  closed  form  analytical  approach  taken  in  references  11-13  is 
not  applicable . 

The  approach  taken  has  included  the  use  of  a  transient  dynamic 
finite  element  code,  modified  for  composites  analysis.  The  code 
selected  for  this  purpose  was  SAP  IV  (ref.  14) .  The  modifications 
made  within  the  finite  element  code  have  allowed  for  the  computation 
of  composite  laminate  properties,  prediction  of  layer  and  interlami¬ 
nar  stresses ,  failure  analysis  and  incorporation  of  predicted  damage 
into  the  subsequent  dynamic  response. 

The  finite  element  method  is  well  suited  for  this  type  of  analy¬ 
sis  since  special  variations  in  material  properties  are  easily  incor¬ 
porated.  Thus,  local  damage  can  be  included  without  affecting  the 
stiffness  of  adjacent  material. 

The  computerized  analysis  procedure,  CLIP  (Composite  Laminate 
Impact  Program),  has  been  utilized  to  predict  the  dynamic  response  of 
a  laminated  plate  subjected  to  low  velocity  impact.  The  predictions 
included  damage  initiation  and  growth  during  time  of  the  dynamic  re¬ 
sponse.  A  description  of  the  code  and  users  guide  are  included  as 
an  appendix  to  this  report. 


II.  ANALYTICAL  METHODOLOGY 


The  analysis  of  laminated  composites  subjected  to  impact  load¬ 
ings  consists  of  several  distinct  procedures.  The  dynamic  response  v 

of  the  system  is  required.  Coupled  with  this  are  laminate  stress 
and  failure  analyses.  Additionally,  predicted  damage  must  be  incor¬ 
porated  back  into  the  dynamic  response  predictions.  The  analysis 
operates  in  a  fashion  where  each  procedure  is  dependent  on  the  others. 

The  dynamic  response  affects  the  stresses.  The  stresses  affect  the 
failure  modes  and  locations.  The  failure  modes  and  locations  in 
turn  affect  the  dynamic  response. 

Each  of  these  areas  and  their  effects  upon  the  other  procedures 
are  discussed  here. 

DYNAMIC  ANALYSIS 

The  dynamic  analysis  routines  used  in  developing  the  CLIP  code 
were  taken  from  SAP  IV.  The  routines  v/hich  were  utilized  include  an 
anisotropic  thin  shell  finite  element,  capable  of  bending  and  mem¬ 
brane  loading,  and  a  time  integration  scheme  which  integrates  the 
equations  of  motion  to  predict  the  time  dependent  structural  response. 

Thin  Shell  Finite  Element 

The  finite  element  used  in  the  analysis  combines  a  bending 
element  with  a  linear  curvature  field  and  a  membrane  element  with 
a  constant  strain  field.  The  two  elements  are  combined  in  such  a 
fashion  that  bending-extensional  coupling  cannot  be  modeled.  Thus, 
it  is  required  that  laminates  modeled  be  mid-plane  symmetric. 

The  element  as  formulated  in  SAP  IV  used  the  same  material  prop¬ 
erties  for  bending  and  extension.  This  has  been  modified  such  that 

« 

nonhomogeneous  materials  can  be  properly  modeled  with  different  prop¬ 
erties  in  extension  and  bending.  The  shell  element  does  have  the 
capability  to  model  both  extensional-shear  coupling  and  bending- 


4 


twisting  coupling  as  it  utilizes  a  full  plane  stress  stiffness  ma¬ 
trix  in  its  foimiulation  and,  as  such,  is  well  suited  for  composites 
analysis. 

For  dynamic  analysis,  it  is  required  that  a  mass  matrix  also  be 
formulated.  The  method  used  in  formulating  the  thin  shell  element 
consists  of  generating  a  lumped  mass  vector.  Hence,  only  diagonal 
mass  terms  are  formed.  Additionally,  no  rotary  inertial  terms  are 
included.  Thus,  the  mass  is  distributed  at  translationary  degrees 
of  freedom  only. 

Since  dynamic  responses  typically  may  include  the  effects  of 
material  damping,  some  provision  must  be  made  for  these  effects  in 
an  analysis.  The  method  used  in  SAP  IV  consists  of  Raleigh  damping. 
This  is  a  convenient  formulation  as  it  does  not  require  the  formula¬ 
tion  of  elemental  damping  matrices.  The  effects  of  damping  are  in¬ 
cluded  at  the  global  mass  and  stiffness  level  and,  hence,  will  be 
discussed  in  the  next  section. 

Time  Integration  Routines 

The  transient  time  analysis  is  performed  using  the  Wilson-0 
method  (refs.  14,  15).  This  procedure  is  a  modification  of  the  lin¬ 
ear  acceleration  scheme.  The  method  is  unconditionally  numerically 
stable  for  any  choice  of  time  step.  The  results  of  the  analysis  are 
dependent  on  the  time  step,  however,  due  to  numerical  damping. 

The  effects  of  numerical  damping  can  be  overcome  by  judicious 
time  step  selection.  The  approach  involves  determining  the  highest 
structural  frequency  of  interest  and  then  selecting  a  time  step  which 
is  a  small  fraction  of  the  period  of  this  frequency.  In  reference 
15,  it  is  shown  that  for  the  Wilson-0  method  an  amplitude  decay  of 
one  percent  per  cycle  can  be  expected  for  a  time  step  to  bending 
mode  period  ratio  of  approximately  0.045. 

The  scheme  used  for  material  damping  in  the  analysis  consists 
of  Raleigh  damping.  In  this  method,  the  damping  is  assumed  to  be 
proportional  to  both  the  stiffness  and  mass  matrices.  This  is  con¬ 
venient  since  no  elemental  damping  matrices  are  required.  Also, 


since  typically  the  structure  to  be  analyzed  with  CLIP  will  be  of 
one  material,  this  type  of  damping  is  most  suitable. 

Initial  Conditions 


The  impact  analysis  performed  by  CLIP  assumes  a  perfectly  plas¬ 
tic  impact.  This  implies  that  the  impacting  mass  attaches  to  the 
plate  and  remains  in  contact.  This  type  of  impact  also  implies  a 
conservation  of  momentum  within  the  system. 

These  effects  are  incorporated  by  lumping  the  impact  mass  at 
specified  nodes  on  the  plate  structure.  An  initial  velocity  is  then 
applied  at  these  nodes.  The  velocity  used  is  scaled  such  that  the 
product  of  impact  mass  and  impactor  velocity  is  equal  to  the  product 
of  impact  mass  plus  nodal  masses  and  the  applied  initial  velocity, 
thereby  conserving  the  impact  momentum. 

There  are  provisions  which  have  been  added  to  allow  for  the  im¬ 
pacting  mass  to  rebound  from  the  plate  structure.  The  contact  force 
between  the  impactor  and  plate  is  computed  based  on  the  product  of 
the  accelerations  of  the  impacted  nodes  and  the  impactor  mass.  When 
this  force  becomes  tensile,  the  mass  detaches  and  the  analysis  be¬ 
comes  a  free  vibration  problem. 

COMPOSITE  AND  STRESS  ANALYSIS 

The  purpose  of  the  composites  analysis  is  two-fold.  First,  lami¬ 
nate  properties  must  be  generated  for  input  to  the  thin  shell  finite 
element  and  second,  composite  inter-  and  intralaminar  stresses  must 
be  predicted  for  use  with  a  failure  analysis.  The  generation  of 
properties  and  in-plane  stress  analysis  of  undamaged  finite  elements 
directly  follows  classical  lamination  theory  and  as  such  will  not  be 
detailed  here.  The  areas  which  need  explanation  include  the  genera¬ 
tion  of  stress,  moment,  and  transverse  shear  resultants  from  the  nodal 
forces  and  moments,  and  the  prediction  of  transverse  (interlaminar) 
shear  stresses  from  the  shear  resultants . 
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The  stress  recovery  portions  of  the  thin  shell  finite  element 
in  SAP  IV  are  designed  to  generate  both  stresses  and  moment  resul¬ 
tants  within  a  triangular  or  arbitrary  quadrilateral  element.  The 
capability  for  predicting  transverse  shear  resultants  did  not 
exist,  however.  In  order  to  include  this  capability  and  minimize 
redundancy  within  the  CLIP  code,  the  SAP  IV  stress  recovery  routines 
were  removed.  The  shell  element  procedures  were  then  modified  to 
produce  nodal  forces  and  moments  directly  from  the  elemental  stiff¬ 
ness  matrices.  By  using  equilibrium  considerations,  stress,  moment 
and  transverse  shear  resultants  are  generated  directly  from  these 
nodal  forces  and  moments.  These  considerations  are  detailed  in  Ap¬ 
pendix  A. 

As  a  consequence  of  the  required  removal  of  the  SAP  IV  stress 
recovery  routines,  and  the  inclusion  of  the  alternate  method  used, 
all  shell  elements  within  the  CLIP  code  are  required  to  be  rectangu¬ 
lar.  Additionally,  all  elements  must  be  aligned  with  the  global 
coordinate  system. 

It  is  Still  possible  to  use  nonrectangular  quadrilateral  ele¬ 
ments  and  obtain  the  correct  displacement  response  but  the  subsequent 
stress  analysis  would  be  in  error. 

The  prediction  of  transverse  (interlaminar)  shear  stress  follows 
closely  the  method  used  in  classical  strength  of  materials  texts  for 
beam  bending  shear  stresses.  The  finite  element  models  used  do  not 
include  shear  deformation  and,  therefore,  the  transverse  shear  re¬ 
sultants  are  determined  through  equilibrium.  This  method  produces 
satisfactory  results  if  the  thickness  of  the  plate  is  small  and  the 
bending  response  dominates  the  overall  deflections. 

The  method  used  for  interlaminar  shear  stress  calculation  in 
both  intact  and  damaged  elements  is  detailed  in  Appendix  B.  Pre¬ 
dicting  interlaminar  shear  stresses  in  damaged  elements  is  consider¬ 
ably  more  difficult  and  will  be  discussed  in  a  later  section. 

The  calculation  of  interlaminar  normal  stresses  within  the  lami¬ 
nated  plate  is  possible  through  the  use  of  stress  equilibrium.  How¬ 
ever,  these  stresses  should  have  significant  magnitude  only  directly 
under  the  impacting  mass.  Thus,  the  effort  required  to  do  a  rigorous 
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analysis  was  deemed  too  extreme.  As  a  first  approximation,  inter¬ 
laminar  normal  stresses  are  computed  directly  under  the  impact  site 
as  a  linear  function  of  the  thickness.  The  contact  force  is  dis¬ 
tributed  over  the  affected  elements  and  scaled  such  that  the  stress 
on  the  back  face  is  zero . 

Failure  Analysis 

As  the  inter-  and  intralaminar  composite  stresses  are  computed, 
it  is  necessary  to  evaluate  whether  or  not  failure  has  occurred. 

This  evaluation  must  predict  both  the  presence  of  failure  and  the 
type  or  mode  of  failure.  The  mode  of  failure  is  required  such  that 
the  specific  elastic  properties  affected  can  be  modified  without  af¬ 
fecting  the  other  material  parameters. 

The  failure  criteria  selected  for  use  in  the  CLIP  code  are  listed 
in  tcible  1.  These  criteria  are  applied  ply  by  ply  for  in-plane  fail¬ 
ure  analysis,  and  interface  by  interface  for  interlaminar  failure 
analysis.  The  components  of  stress  used  in  the  failure  analysis  are 
described  in  figure  1 . 

Incorporation  of  Predicted  Damage 

The  primary  goal  of  the  current  study  is  to  track  damage  accvimu- 
lation  and  growth  throughout  the  impact  event.  This  involves  calcu¬ 
lating  damage  at  selected  time  steps  and  incorporating  the  effects  of 
this  damage  back  into  the  dynamic  solution.  The  failure  criteria  de¬ 
scribed  in  the  previous  section  are  used  to  calculate  the  modes  of 
failure  within  the  shell  finite  elements.  Wlien  the  modes  of  damage 
are  found  in  an  element,  its  stiffness  matrix  is  reformed  and  the 
new  element  is  incorporated  into  the  dynamic  analysis.  This  proce¬ 
dure  is  accomplished  by  subtracting  out  the  original  element  from 
the  unfactored  global  stiffness  matrix  and  adding  the  new,  damaged 
element.  This,  of  course,  requires  that  the  global  stiffness  matrix 
be  decomposed  again. 
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Ply  Damage 


When  the  damage  mode  predicted  consists  of  ply  damage,  it  is 
necessary  to  compute  new  laminate  stiffnesses  for  bending,  extension 
and  bending-extensional  coupling.  The  new  bending  and  extensional 
properties  are  input  back  into  the  shell  finite  element  routines  for 
formulating  the  new  elemental  stiffness  matrix.  The  bending- 
extensional  coupling  terms  are  used  only  in  the  stress  recovery  in¬ 
formation  since,  as  was  mentioned  before,  this  type  of  material  be¬ 
havior  cannot  be  modeled  in  the  displacement  solution. 

The  way  in  which  the  laminate  properties  are  changed  depends 
upon  the  type  of  ply  damage.  If  the  ply  damage  includes  fiber  fail¬ 
ure,  then  the  entire  ply  is  removed  from  the  laminate  model.  If, 
however,  matrix  failure  is  the  only  failure  mode,  it  is  assumed  that 
the  fiber  can  still  carry  load.  Thus,  only  the  transverse  properties 
of  the  ply  are  removed.  These  laminate  modifications  apply  to  the 
specific  element  under  consideration-  Thus,  if  only  one  element 
sustains  damage,  the  other  elements  are  not  affected. 

When  an  interior  ply  fails  under  the  dynamic  loading,  it  is 
assumed  that  the  constraint  of  the  adjacent  material  will  cause  the 
strain  distribution  to  remain  linear  through  the  laminate  thickness. 
Because  of  this,  it  is  not  necessary  to  compute  properties  for  two 
or  more  separate  laminates  for  the  damaged  shell  elements.  This 
type  of  damage  does,  however,  pose  some  difficulties  in  finding  the 
interlaminar  stresses.  The  methodology  used  for  predicting  inter¬ 
laminar  shear  stresses  is  described  in  Appendix  B. 

Within  the  CLIP  code,  the  types  and  locations  of  damage  are 
saved  for  each  damaged  element.  This  information  is  then  used  to 
evaluate  subsequent  damage  modes.  Under  certain  circumstances,  it 
might  be  possible  for  the  numerical  procedure  to  predict  the  same 
damage  mode  and  location  more  than  once.  This  is  physically  unreal¬ 
istic  though,  and  if  the  CLIP  code  senses  this,  the  analysis  is 
stopped. 
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Interlaminar  Damage 

When  interlaminar  damage  is  predicted  due  to  the  dynamic  load¬ 
ing,  the  procedure  used  is  different  from  that  for  ply  damage.  Be¬ 
cause  of  the  assumption  that  adjacent  material  restrains  the  curv^ 
ature  and  strain  distribution  when  local  damage  is  present,  the  only 
effect  of  a  delamination  is  an  adjustirient  to  the  transverse  shear 
stress  field.  The  shell  finite  element  does  not  consider  shear 
deformation  and,  therefore,  the  increased  shear  deformation  asso¬ 
ciated  with  a  crack  tip  singularity  cannot  be  modeled.  The 
introduction  of  interlaminar  delaminations  does  not  require  a  re¬ 
formulation  of  the  elemental  or  global  stiffness  matrices.  A  full 
description  of  the  assumptions  implicit  in  this  method,  the  effects 
on  the  shear  stress  distribution,  and  the  formulation  of  the  shear 
stress  calculations  is  contained  in  Appendix  B. 

SUMMARY  OF  ANALYSIS  METHODOLOGY 

The  analysis  methodology  described  previously  including  Appen¬ 
dices  A  and  B  detail  the  approach  taken  for  analyzing  low  velocity 
impact  of  composite  plates.  The  assimiptions  implicit  in  this  metho¬ 
dology,  as  well  as  the  capabilities  and  limitations  of  the  analysis, 
are  stunmarized  here  for  clarity. 

Analysis  Assumptions 
1.  Small  Deflection  Analysis 

The  displacements  of  the  plate  structure  are  sufficiently 
small  such  that  the  original  geometry  is  applicable  throughout 
the  analysis. 
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2.  Displacement  Fields 


Only  bending  and  membrane  displacements  are  modeled. 

The  effects  of  shear  deformation  are  insignificant  with  re¬ 
spect  to  bending  deformations.  Typically,  this  condition  is 
satisfied  if  the  plate  thickness  is  small  compared  to  the  other 
plate  dimensions. 

3.  Material  Properties 

Static  material  constants  and  strengths  are  used.  The 
effects  of  time-dependent  material  properties  are  insignifi¬ 
cant  in  a  realistic  structural  laminate.  This  is  primarily 
due  to  the  presence  of  fibers  in  many  directions. 

4.  Superposition  of  Static  and  Dynamic  Displacements 

Static  pre-stress  displacement  fields  are  added  directly 
to  the  dynamic  displacements.  No  stability  analysis  is  per¬ 
formed  as  a  consequence  of  (1) . 

5.  Laminate  Modeling 

Bending-extensional  coupling  is  not  included  in  the 
analysis.  This  effectively  requires  that  all  laminates 
modeled  be  mid-plane  symmetric.  The  assumption  is  made  that 
damage  induced  bending-extensional  coupling  will  be  highly 
localized  and,  therefore,  negligible. 

6.  Perfectly  Plastic  Impact 

The  impact  involves  a  complete  momentum  transfer  from 
the  impacting  mass  to  the  composite  plate  structure.  This 


also  implies  that  contact  effects  (Hertzian  Contact)  and 
wave  propagation  effects  are  negligible,,  Additionally, 
since  contact  effects  are  not  modeled,  the  impact  mass 
can  only  rebound  due  to  plate  dynamics, 

7,  Interlaminar  Normal  Stresses 

Interlaminar  normal  stresses  are  computed  only  di- 
rectly  under  the  impacting  mass,  A  linear  approximation 
is  used  through  the  thickness  of  the  plate  such  that  the 
back  surface  has  zero  stress  and  the  impacted  surface 
carries  the  contact  force  distributed  over  all  affected 
elements. 


8,  Matrix  Ply  Damage 

Damage  produced  by  combined  <^22'  ^12  fields. 

Matrix  dominated  ply  moduli  E2  and  G^2  zero 

in  the  affected  ply  within  the  affected  finite  element. 


9.  Fiber  Ply  Damage 

Damage  produced  by  stresses .  All  lamina  moduli 

are  set  to  zero  in  the  affected  ply  within  the  affected 
finite  element. 

10.  Interlaminar  Delamination 

No  stiffness  effects  as  a  direct  consequence  of  (2) . 
Increased  shear  deformations  associated  with  a  singular 
shear  stress  distribution  at  the  crack  tip  are  not  modeled. 


Damage  Stress  Effects 


11.  Matrix  Ply  Damage 

Lamina  stresses  022  automatically  set  to  zero 

as  a  consequence  of  (8) .  Additionally,  interlaminar  stresses 

a  and  a  are  set  to  zero  within  the  damaged  ply. 
xz  yz  ^  f  1 

12.  Fiber  Ply  Damage 


Lamina  stresses  cr^^,  022  <^2.2  automatically  set  to 

zero  as  a  consequence  of  (S) .  Additionally,  interlaminar 
stresses  and  are  set  to  zero  within  the  damaged  ply. 


13.  Interlaminar  Delamination 


Interlaminar  stress  components  and  are  set  to 

zero  at  the  affected  ply  interface  within  the  affected  finite 
element. 


Damage  Propagation 


14.  Ply  Damage 

Damage  predicted  in  plies  propagated  due  to  load  redis¬ 
tribution.  The  load  distribution  is  changed  because  of  local 
stiffness  changes  (8) ,  (9) . 

15.  Interlaminar  Delamination 

Delaminations  do  not  propagate  to  adjacent  elements.  De¬ 
laminations  may  occur  in  adjacent  elements  due  to  shear  force 
distributions  but  these  distributions  are  not  changed  as  a 
result  of  interlaminar  damage  (10)  . 
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III.  LOW  VELOCITY  IMPACT  ANALYSIS 


The  analysis  methodology  described  in  the  previous  sections  has 
been  utilized  to  predict  the  response  of  a  clamped  rectangular  lami¬ 
nated  composite  plate  subjected  to  a  low  velocity,  low  mass  impact. 
The  analyses  performed  have  provided  information  relating  to  the 
modes  and  locations  of  impact  induced  damage,  the  displacement  re¬ 
sponse  of  the  composite  plate  and  the  contact  force  between  the  im.-" 
pacting  mass  and  the  plate  structure.  This  information  is  generated 
at  various  times  throughout  the  duration  of  the  impact  event. 

In  the  course  of  predicting  the  response  of  the  plate,  it  was 
required  to  make  several  computer  analyses..  These  included  solutions 
with  the  entire  impact  mass  lumped  at  a  single,  central  node  and  so¬ 
lutions  with  the  impacting  mass  distributed  over  a  group  of  centrally 
located  nodes. 

The  laminate  configuration,  impact  parameters,  finite  element 
models  and  the  results  of  the  various  analyses  are  described  and 
discussed  here. 

ANALYSIS  PARAMETERS 

Before  discussing  the  various  analyses  performed,  the  various 
parameters  relating  to  the  materials ,  finite  element  models  and  im¬ 
pact  conditions  need  to  be  described.  The  analyses  performed  were 
part  of  an  effort  to  model  one  of  many  impact  experiments  which  have 
been  performed  at  NASA  Langley.  Hence,  the  materials  and  configura¬ 
tions  correspond  to  this  experiment . 

Laminate  Configuration 


The  laminate  selected  for  the  impact  analysis  was  an  eight  ply, 
quasi-isotropic  configuration.  The  stacking  sequence  analyzed  v/as 
t45/0/-45/90] s •  The  material  properties  used  correspond  to  a  T300/ 
5208  Gr/Ep  system  and  are  listed  in  table  2.  The  unidirectional 
properties  used  are  typical  static  data  and  were  taken  from  refer- 
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ence  16.  It  should  be  noted  that  variation  in  these  properties 
could  greatly  affect  the  solution.  Therefore,  a  design  application 
of  the  analysis  would  require  the  characterization  of  the  material 
under  consideration  to  avoid  the  wide  variation  in  material  data  re¬ 
ported  in  the  literature. 

The  ply  thickness  used  was  0.01321  cm,  yielding  a  total  laminate 
thickness  of  0.1056  cm.  This  laminate  configuration  exhibits  bend¬ 
ing-twisting  coupling  which  must  be  taken  into  account  in  the  finite 
element  model. 


Finite  Element  Models 

The  finite  element  model  used  for  most  of  the  analyses  made  is 
shown  in  figure  2.  The  model  encompasses  the  entire  plate  structure 
as  required  by  the  bending-twisting  coupling  present  in  the  laminate 
to  be  analyzed.  The  model  dimensions  are  10.16  cm.  by  15,24  cm. 

In  figure  2,  the  shaded  area  represents  the  elements  which  were 
selected  for  stress  analysis.  The  CLIP  code  has  been  developed  such 
that  only  specified  elements  have  stress  calculations  performed  (see 
Appendix  C) .  This  was  done  in  an  effort  to  maximize  computational 
economy . 

The  central  section  of  the  plate  was  chosen  for  stress  calcula¬ 
tions  as  this  is  the  area  where  any  impact  induced  damage  should  oc¬ 
cur,  More  elements  could  have  been  selected  for  stress  analysis  but 
this  was  deemed  unnecessary  since  the  impact  site  was  directly  at 
the  center  of  the  plate. 

The  model  was  developed  to  represent  a  clamped  plate  and  as 
such,  the  edge  nodes  are  constrained  against  all  rotations  and  dis¬ 
placements.  The  model  contains  704  elements,  759  nodes  and  1953 
active  degrees  of  freedom.  The  number  of  active  degrees  of  freedom 
is  minimized  by  constraining  all  in-plane  displacements. 

In  the  finite  element  model,  the  maxim\im  element  aspect  ratio 
is  eight.  Within  the  area  where  stress  calculations  are  performed, 
the  maximiira  element  aspect  ratio  is  4 , 
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In  order  to  verify  the  validity  of  this  model,  a  static  solu¬ 
tion  was  run.  The  loading  consisted  of  a  uniform  pressure  distrib¬ 
uted  over  the  entire  plate  surface.  The  results  of  this  solution 
were  then  compared  with  a  one-term  Ritz  solution  taken  from  refer¬ 
ence  17.  For  this  comparison,  it  was  necessary  to  assume  that  the 
quasi-isotropic  [45/0/-45/90] g  laminate  behaves  as  a  specially  or¬ 
thotropic  laminate.  Hence,  the  bending- twisting  coupling  terms, 

Di6f  ^26  ignored.  Even  with  this  approximation,  the  two  solu¬ 

tions  compare  within  3.3%.  Hence,  the  model  has  been  shown  to  be 
valid. 

During  the  course  of  performing  the  impact  analyses,  it  was 
necessary  to  develop  another  finite  element  model.  This  model  is 
shown  in  figure  3.  The  only  difference  between  the  two  models  is 
the  removal  of  four  elements  and  one  node  at  the  center  of  the  plate 
in  the  second  model.  The  rationale  for  developing  this  second  model 
will  be  discussed  later  in  this  report. 

Impact  Parameters 

The  impact  mass  and  velocity  modeled  in  the  analyses  are  listed 
in  table  3.  All  of  the  analyses  used  these  parameters.  The  vari¬ 
ous  analyses  did  in  some  cases  represent  different  distributions  of 
the  impact  mass  on  the  plate,  however.  The  different  distributions 
used  are  depicted  in  figure  4.  Each  of  the  sections  shown  are  repre¬ 
sentative  of  the  geometric  center  of  the  plate.  Obviously,  the 
third  impact  mass  distribution  shown  in  figure  4  was  utilized  with 
the  second  finite  element  mesh  which  had  the  four  central  elements 
removed . 


Integration  Time  Step 

The  selection  of  the  integration  time  step  is  one  of  the  more 
critical  steps  in  defining  the  impact  model.  The  time  step  selected 
must  be  small  enough  to  adequately  determine  the  response  of  all 
critical  bending  modes  while  not  requiring  an  excessive  number  of 
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time  steps  to  investigate  the  impact  event.  The  approach  is  to  com¬ 
pute  the  periods  of  the  plate  natural  frequencies  and  determine  which 
are  critical. 

An  analysis  was  made  using  a  solution  given  in  reference  17  to 
determine  the  natural  frequencies  for  the  clamped  plate  in  question. 
As  before,  it  was  necessary  to  model  the  plate  without  the  bending¬ 
twisting  coupling  terms.  This  implies  that  all  frequencies  calcula¬ 
ted  are  higher  than  actual  since  neglecting  the  coupling  terms  ef¬ 
fectively  increases  the  plate  bending  stiffnesses.  Additionally, 
the  predictions  were  made  without  including  the  effects  of  impact 
induced  damage.  Thus,  the  stiffnesses  were  again  higher  in  the 
predictions  than  can  be  expected  in  the  impact  analyses.  Hence,  the 
computed  frequencies  can  be  expected  to  be  somewhat  higher  than  the 
actual  natural  frequencies 

Based  on  these  calculations,  a  time  step  of  1  ysec  was  selected 
for  the  initial  impact  analysis.  This  corresponds  to  approximately 
1/19  of  the  10-10  bending  mode  period.  The  results  of  the  initial 
impact  analysis  indicate  that  the  1  ysec  time  step  was  considerably 
smaller  than  required.  The  bending  shapes  of  the  plate  did  not  in¬ 
volve  frequencies  this  large  and,  hence,  for  the  remaining  analyses, 
a  time  step  of  2.5  ysec  was  used  with  no  apparent  degradation  of  the 
results . 

In  conjunction  with  the  time  step  selection,  it  must  be  deter¬ 
mined  how  often  to  perform  stress  calculations.  In  each  of  the 
analyses  made  where  stress  calculations  were  included,  the  frequency 
of  stress  calculation  was  every  five  time  steps.  In  terms  of  com¬ 
puting  damage  growth,  the  optimum  frequency  of  stress  calculation 
would  be  to  compute  them  every  time  step.  Considerations  of  the 
cost  and  time  involved  precluded  this,  however. 


IMPACT  ANALYSES 

In  order  to  predict  the.  response  of  the  laminated  plate  described 
when  subjected  to  the  impact  conditions  also  described,  it  was  neces¬ 
sary  to  perform  four  separate  analyses.  The  first  three  analyses 
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were  performed  with  the  finite  element  model  shown  in  figure  2. 

The  last  analysis  was  made  using  the  model  in  figure  3. 

The  first  analysis  was  made  without  the  inclusion  of  stress 
calculations.  This  analysis  was  performed  in  order  to  verify 
the  time  step  selection  and  determine  the  duration  of  the  impact 
event.  (See  Appendix  C.) 

The  second  and  third  analyses  were  full  analyses  including 
stress  analysis.  In  the  second  analysis,  the  impact  mass  was  lumped 
at  the  center  node  of  the  finite  element  model.  The  results  of  this 
analysis  prompted  the  distributed  impact  mass  utilized  in  the  third 
analysis.  In  both  of  these  solutions,  the  damage  computed  was  so 
extensive  that  computerized  procedure  terminated  the  solution. 

The  fourth  solution  was  made  in  an  attempt  to  model  the  dam¬ 
age  growth  beyond  the  point  at  which  the  computer  program  had 
terminated  in  the  second  and  third  solutions.  Hence,  the  removal 
of  the  four  central  elements  in  figure  3. 

Convenient  groupings  of  the  damage  predictions  of  these  analy¬ 
ses  can  be  found  in  table  4.  The  information  in  table  4  will  aid 
in  comparisons  of  the  growth  of  damage  between  solutions  and  in 
comparisons  of  the  types  of  damage  within  each  solution. 

Displacement  Solution 

The  first  analysis  performed  was  simply  a  dynamic  displacement 
response  solution.  The  impact  mass  was  lumped  at  the  center  plate 
node.  The  integration  time  step  was  1,0  ysec. 

The  displacement  response  of  the  laminated  plate  plotted 
through  the  center  of  the  plate,  along  the  axis  corresponding  to 
the  smaller  dimension  of  the  plate,  is  shown  in  figures  5  and  6. 

The  displacements  in  figure  5  represent  the  very  short  time  response 
of  the  plate  while  figure  6  depicts  the  longer  time  response. 

Comparing  the  two  figures  demonstrates  a  fundamental  difference 
between  the  early  displacement  fields  (fig.  5)  and  later  displace¬ 
ment  fields  (fig.  6) .  At  very  early  times  in  the  impact  event, 
the  displacements  can  be  characterized  as  a  local  phenomenon.  At 
T  =2.5  X  lO”^  seconds,  the  major  displacement  response  can  be  seen 
to  exist  at  the  center  of  the  plate,  with  the  rest  of  the  structure 
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remaining  nearly  motionless.  As  time  progresses,  the  extent  of  the 
plate  with  significant  displacements  can  be  seen  to  be  progressing 
outwards  to  the  edges  of  the  plate. 

The  progress  of  outward  spreading  of  the  major  displacement 
response  is  complete  atT=2.0xlO^  seconds,  as  can  be  seen  in 
figure  6.  All  of  the  displacement  fields  in  figure  6  can  be  char¬ 
acterized  as  a  predominance  of  the  third  bending  mode  shape.  The 
third  mode  rather  than  the  first  mode  is  excited  due  to  the  clamped 
boundary  conditions . 

Another  aspect  of  the  mode  shapes  excited  relates  to  the  inte¬ 
gration  time  step  selected.  Previously  it  was  stated  that  the 
1  ysec.  time  step  was  sufficient  to  characterize  the  10-10  bending 
mode  shape  of  the  impacted  plate.  It  is  quite  obvious  from  figure 
5  that  the  mode  shapes  present  do  not  approach  the  10-10  mode. 

Hence,  for  all  remaining  analyses,  the  time  step  was  increased  to 
2.5  ysec. 

In  setting  up  the  analysis,  two  unfortunate  situations  were  in¬ 
cluded.  First,  the  printing  of  displacements  was  set  up  in  such  a 
fashion  that  it  was  not  possible  to  observe  the  displacement  fields 
along  the  longer  dimension  of  the  plate.  This  is  apparently  of  lit¬ 
tle  consequence  since  the  displacements  of  figures  5  and  6  have  pro¬ 
vided  sufficient  information.  The  second  problem  involves  the 
amount  of  time  required  for  the  duration  of  the  impact  event.  The 
solution  was  set  up  with  five  hundred  time  steps.  Hence,  the  total 
time  allowed  was  5  x  10“4  seconds.  At  the  end  of  this  time,  and 
therefore  the  end  of  the  solution,  the  impact  mass  had  not  rebounded 
and  the  maximum  displacements  had  not  yet  been  reached.  This  caused 
some  difficulty  since  part  of  the  rationale  for  performing  this  solu 
tion  involved  determining  the  time  duration  of  the  impact  event. 

In  order  to  make  an  estimate  of  the  duration  of  the  impact,  it 
was  necessary  to  consider  the  contact  force  between  the  impact  mass 
and  the  plate.  In  figure  7,  a  plot  of  the  force-time  response  of 
the  impact  event  is  shown.  The  force  appears  to  be  somewhat  erratic 
as  it  is  computed  as  the  product  of  the  impact  mass  and  the  accelera 
tion  of  the  plate  node  where  the  mass  is  lumped.  The  accelerations 
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of  the  node  in  question  reflect  considerably  higher  frequencies  than 
the  displacements  and,  hence,  the  contact  force  calculated  also  con¬ 
tains  the  high  bending  mode  frequencies. 

Evaluations  of  the  minimum  value  of  the  forces  shown  in  figure 

7  indicate  that  the  impact  mass  was  close  to  rebounding  when  the 

analysis  terminated.  Additionally,  a  displacement  analysis  made 

with  a  less  refined  finite  element  model  indicated  that  the  maximum 

displacement  value  for  these  impact  conditions  was  nearly  attained 

at  0.5  msec.  Therefore,  for  the  remainder  of  the  impact  analyses, 

the  maximum  time  allowed  for  the  solution  was  increased  to  8.75  x 
-4 

10  seconds.  This  corresponds  to  350  time  steps  at  the  new  inte¬ 
gration  time  step  of  2.5  ysec. 

Stress  Solution,  Mass  Lumped  at  One  Node 

Using  the  new  time  step  described,  a  full  impact  analysis  with 
stress  solution  was  performed.  This  solution  progressed  until  two 
elemental  laminate  stiffness  matrices  became  singular  due  to  damage 
in  all  plies  at  the  end  of  thirty  time  steps. 

The  displacement  time  response  through  the  center  of  the  plate 
along  the  shorter  dimension  of  the  plate  is  shown  in  figure  8.  Com¬ 
paring  these  displacements  with  those  shown  in  figure  5  demonstrates 
little,  if  any,  difference  between  the  two  solutions.  Since  in 
figure  8  considerable  damage  is  present,  especially  at  7.5  x  10 
seconds,  one  would  expect  significant  differences  to  be  present. 

The  reason  that  these  differences  are  not  present  is  simply  that 
the  response  shown  in  both  figures  5  and  8  is  primarily  inertial. 

Not  enough  time  has  passed  for  the  effects  of  the  induced  damage  to 
significantly  affect  the  solution.  Had  the  stress  solutions  contin¬ 
ued  for  a  longer  time,  the  differences  would  have  become  significant. 

In  setting  up  this  solution,  sufficient  displacement  printing 
was  specified  such  that  the  response  of  the  plate  along  the  longer 
axis  could  be  observed.  The  displacement  time  response  through  the 
center  of  the  plate  along  the  longer  dimension  of  the  plate  is  shown 


20 


in  figure  9 .  The  displacement  fields  along  the  long  axis  of  the 
plate  show  considerable  similarity  to  those  along  the  shorter  axis . 

The  primary  difference  is  that  the  changes  in  slope  of  the  curves 
are  slightly  more  gradual  along  the  longer  axis.  This  was  expected 
since  the  plate  is  longer  and,  hence,  more  flexible  in  this  direction. 

In  both  figures  8  and  9,  the  most  striking  result  is  the  limited 
area  of  the  plate  with  significant  displacement  response  during  the 
early  moments  of  the  impact  event.  As  time  progresses,  the  displace¬ 
ment  "waves"  move  outward  until  the  entire  plate  is  affected  and  one 
would  expect  high  moments  and  shears  corresponding  to  the  compressed 
displacement  fields . 

The  moment  resultants  through  the  plate  center  along  the  shorter 
and  longer  plate  axes  are  shown  in  figures  10  and  11  respectively ^ 

The  moment  resultants  correspond  to  T  =  5.0  x  10~5  seconds  and, 
hence,  the  middle  displacement  fields  in  figures  8  and  9.  The  moment 
resultants  depicted  in  figures  10  and  11  represent  averages  of  the 
elements  on  either  side  of  the  plate  centerlines.  The  specific 
values  vary  slightly  on  either  side  of  the  centerlines  due  to  the 
bending/ twisting  coupling  described  previously.  This  coupling  also 
produces  twisting  moment  resultants.  The  twisting  moments  are  not 
shown,  however,  since  they  are  of  such  small  magnitude.  The  maxi¬ 
mum  twisting  moment  predicted  was  less  than  8%  of  the  peak  value. 

The  moment  resultants  shown  in  figures  10  and  11  demonstrate 
that  the  major  effect  of  the  impact  is  highly  localized  at  5.0  x  10“^ 
seconds.  This  is  in  agreement  with  the  displacement  fields  depicted 
in  figures  8  and  9 . 

The  transverse  shear  resultants  corresponding  to  the  moment 
resultants  are  shown  in  figures  12  and  13.  Once  again,  the  major 
loading  occurs  in  a  local  region  surrounding  the  impact  site.  This 
is  due,  of  course,  to  the  large  moment  resultant  gradients  at  the 
plate  center  (figs.  10  and  11) . 

The  transverse  shear  resultants  are  plotted  in  the  same  fashion 
as  the  moment  resultants.  The  figures  represent  averages  of  the 
two  elements  on  either  side  of  the  plate  centerlines.  The  shear 
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resultants  in  individual  elements  also  contain  the  effects  of  the 
bending- twisting  coupling  and,  therefore,  are  slightly  different  on 
either  side  of  the  plate  centerlines. 

The  moment  and  transverse  shear  resultants  depicted  in  figures 
10-13  are  sufficiently  large  to  produce  significant  damage  at  the 
center  of  the  plate.  At  the  end  of  five  integration  steps,  a  stress 
analysis  was  performed.  The  elements  which  suffered  in-plane  and 
interlaminar  damage  at  this  time  are  shown  in  figures  14  and  15 
respectively.  The  grid  on  which  the  damage  is  shown  encompasses 
all  elements  which  were  selected  for  stress  analysis. 

The  ply  damage  shown  in  figure  14  consists  primarily  of  trans¬ 
verse  (matrix)  cracking  in  the  bottom  45“  ply.  Six  of  the  elements 
also  experience  matrix  cracking  within  the  next  interior  ply  (0“) . 
Thus,  the  region  shown  in  figure  14  has  sustained  considerable  damage 
in  terms  of  surface  area  but  little  through  the  thickness. 

The  interlaminar  damage  shown  in  figure  15  is  more  interesting. 
Each  of  the  four  elements  depicted  has  experienced  delamination  at 
the  five  innermost  interfaces.  This  damage  is  quite  extensive 
and  indicates  that  the  transverse  shear  resultants  are  extremely 
high  in  this  region  and  at  this  time.  In  practical  terms,  the 
material  in  this  region  must  be  considered  fully  degraded. 

The  most  interesting  feature  of  the  interlaminar  damage  is 
its  presence.  Had  one  simply  applied  a  static  point  load  at  the 
center  of  the  plate,  no  delamination  would  have  occurred.  However, 
at  the  very  early  time  at  which  the  stresses  were  computed,  extreme¬ 
ly  sharp  bending  gradients  are  present.  This  is  due  to  the  highly 
localized  displacement  response  at  this  early  time,  as  was  shown  at 
later  times  in  figures  8  and  9. 

In  figure  16,  the  elements  with  accumulated  damage  after  2.5 
X  lO”^  seconds  are  depicted.  Comparing  figures  14  and  16,  the  re¬ 
gion  of  damaged  elements  can  be  seen  to  be  growing.  The  majority 
of  the  damaged  elements  in  figure  16  have  one  or  two  back  face  plies 
damaged.  Two  of  the  elements  shown  have  suffered  top  face  damage. 
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These  elements  are  shown  in  figure  17.  Obviously  at  2.5  x  10 
seconds,  the  ply  damage  has  become  extensive,  with  elements  fail¬ 
ing  at  both  outer  surfaces. 

It  is  interesting  that  at  2.5  x  10~^  seconds  no  additional 

interlaminar  damage  is  present.  The  reason  for  this  stems  from 

the  method  used  for  computing  interlaminar  shear  stresses.  These 

stresses  are  computed  at  ply  interfaces  only.  Since  the  remaining 

interface  was  between  0°  and  45“  plies,  the  neutral  surface  will 

not  be  at  the  ply  interface;  the  maximum  shear  stress  will  also 

not  be  at  the  ply  interface.  Thus,  the  analysis  does  not  predict 

the  maximum  shear  stress  in  this  case  and  the  lower  stress  at  the 

interface  may  not  cause  damage. 

At  3.75  X  l0~^  seconds,  the  transverse  shear  resultants  are 

sufficient  to  produce  large  interface  stresses  and  additional  inter 

laminar  damage  is  predicted.  The  new  delamination  is  within  the 
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elements  which  had  sustained  interlaminar  failures  at  1.25  x  10 
seconds.  Thus,  two  of  the  elements  shown  in  figure  15  have  only 
one  remaining  intact  ply  interface.  The  remaining  interface  in 
these  elements  is  between  the  bottom  0°  and  45°  plies. 

In  figures  18,  19  and  20,  the  accumulation  of  ply  damage  dur¬ 
ing  the  remainder  of  the  analysis  is  shown.  No  additional  inter¬ 
laminar  damage  occurred. 

These  figures  show  a  steady  increase  in  the  region  of  impact 
induced  ply  damage.  The  extent  of  the  damage  through  the  thickness 
at  the  end  of  the  analysis  is  depicted  in  figures  20,  21  and  22. 
Figure  20  indicates  the  total  planar  area  of  the  plate  damage.  Fig 
ure  21  shows  which  elements  in  this  region  have  more  than  one  dam¬ 
aged  ply  and  figure  22  shows  elements  which  have  sustained  fiber 
damage. 

The  fiber  breakage  shown  in  figure  22  is  probably  the  most 

critical.  In  each  of  the  elements,  the  first  occurrence  of  fiber 
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damage  was  at  5  x  10  seconds.  Both  elements  sustained  failures 
in  the  bottom  -45°  ply  at  this  time.  At  6.25  x  10  ^  seconds,  ad¬ 
ditional  damage  in  the  bottom  90°  and  0°  plies  occurred.  At  7.5 
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X  10  ^  seconds,  each  ply  in  these  elements  had  sustained  damage  and, 
hence,  the  solution  terminated.  There  was  no  material  left  to  carry 
any  additional  loading. 

The  damage  acciimulated  during  this  solution  was  very  large  in 
its  extent,  both  in  the  planar  dimensions  of  the  plate  and  through 
the  thickness „  The  problem  with  these  results  was  that  the  experi¬ 
ments  carried  out  at  NASA  Langley  did  not  indicate  damage  as  exten¬ 
sive  as  the  computer  analysis  did.  The  experimental  work  generated 
interlaminar  separation  and  transverse  ply  failure  as  the  two  pri¬ 
mary  damage  modes ,  but  apparently  did  not  produce  the  extensive  fiber 
damage.  One  reason  for  this  could  have  been  related  to  the  placement 
of  the  impacting  mass.  Since  the  entire  mass  was  limped  at  one  node 
in  the  analysis,  the  response  of  the  plate  may  have  been  overestimated 
It  was  determined,  therefore,  that  an  additional  analysis  would  be 
made  with  the  impact  mass  distributed  over  five  nodes  rather  than 
limped  at  one  node. 

Stress  Solution,  Mass  Distributed  at  Five  Nodes 

The  impact  mass  distribution  used  in  this  third  analysis  was 
previously  depicted  in  figure  4.  This  distribution  was  chosen  in 
an  effort  to  more  closely  simulate  the  impact  of  a  sphere  on  a  plate. 
The  diameter  encompassed  in  the  distributed  mass  arrangement  corres¬ 
ponds  to  20%  of  the  diameter  of  the  sphere  used  in  the  experiments 
at  NASA  Langley. 

The  displacement  response  predicted  in  this  analysis  was  very 
similar  to  the  previous  analysis.  The  shapes  of  the  curves  were 
nearly  identical  to  those  in  figures  8  and  9.  One  surprising  dif¬ 
ference  was  present,  however.  The  solution  with  the  distributed 
load  produced  larger  peak  displacement  values.  At  2.5  x  10“^  sec¬ 
onds,  the  difference  amounted  to  3%.  When  the  time  increased  to 
5  X  10  ^  seconds,  the  difference  decreased  to  1%. 
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This  result  was  surprising  since  it  was  anticipated  that  the 
spreading  of  the  impact  mass  would  reduce  the  displacements.  The 
reason  for  the  increased  displacements  is  the  initial  velocity  solu¬ 
tion  procedure  used  in  the  analysis.  Since  all  nodes  where  the  im¬ 
pact  mass  is  applied  are  given  an  initial  velocity,  more  of  the  plate 
begins  the  analysis  with  non-zero  velocities.  When  only  one  node  is 
impacted,  the  surrounding  nodes  lag  behind  in  terms  of  velocity. 

When  the  distributed  impact  mass  is  used,  these  surrounding  nodes 
are  excited  at  the  same  velocity.  Thus,  in  the  solution  where  the 
mass  is  lumped  at  one  node,  the  velocity  lag  of  the  surrounding  nodes 
restrains  the  displacements.  This  effect  is  apparently  temporary  in 
the  solution  procedure  as  the  difference  decreases  rapidly  with  time. 

The  distributions  of  damage  in  this  solution  were  somewhat  dif¬ 
ferent  from  those  in  the  previous  solution.  The  extent  of  the  dam¬ 
age  was  nearly  the  same,  however,  and  this  solution  terminated  in  a 
fashion  similar  to  the  previous  case.  The  increased  displacements 
caused  this  solution  to  terminate  earlier  than  the  Ivimped  mass  solu¬ 
tion.  The  procedures  operated  for  20  time  steps  before  the  damage 
was  too  extensive,  while  the  previous  solution  proceeded  for  30  steps. 

At  the  first  stress  calculation,  the  predicted  ply  and  interlami¬ 
nar  damage  is  depicted  in  figures  23  and  24  respectively.  Comparing 
figures  23  and  14,  it  can  be  seen  that  the  extent  of  the  damage  is 
very  similar  for  the  distributed  mass  solution  and  the  lumped  mass 
solution.  The  type  of  damage  is  also  similar  in  that  the  majority 
of  damaged  elements  have  suffered  back  ply  matrix  tensile  failures. 

Comparing  the  interlaminar  damage  predicted  in  the  two  solutions 
(figures  15  and  24)  demonstrates  that  the  distributed  mass  solution 
produces  considerably  more  delamination,  although  the  extent  of  this 
damage  through  the  thickness  of  the  plate  was  comparable  for  the 
two  solutions.  Each  of  the  damaged  elements  in  figure  24  has  suf¬ 
fered  delaminations  between  all  plies  except  the  outer  45° /0°  in¬ 
terfaces.  As  in  the  previous  solution,  this  extensive  delamination 
must  be  considered  as  total  failure  in  these  elements.  The  solu¬ 
tion  procedure  continued,  however,  since  most  of  the  individual  plies 
remained  intact. 
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The  damage  growth  during  the  remainder  of  the  analysis  is  de¬ 
picted  in  figures  25  through  27.  With  each  stress  calculation,  the 
number  of  damaged  elements  grows,  as  does  the  extent  of  the  dam¬ 
age  within  the  elements.  Comparing  figures  27  and  19  indicates 
that  the  number  of  damaged  elements  at  5  x  10  ^  seconds  is  identical 
for  the  two  solutions.  The  types  of  damage  within  the  elements  are 
very  different,  however.  The  damage  in  the  central  four  elements 
(fig.  28)  includes  considerable  fiber  breakage.  The  solution  ter¬ 
minated  because  two  of  these  central  elements  had  sustained  fiber 
damage  in  all  plies  except  the  upper  0°  ply.  This  ply  had  suffered 
matrix  damage,  however.  Thus,  no  material  remained  for  carrying 
the  load. 

The  delamination  predicted  in  the  distributed  mass  solution 
never  progressed  beyond  the  twelve  elements  shown  in  figure  24. 

The  amount  of  delamination  within  these  elements  did  increase,  how- 
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ever.  At  2.5  x  10  seconds,  the  four  central  elements  (fig.  29) 

had  no  remaining  ply  interfaces.  During  the  subsequent  stress  cal- 
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culation  at  3.75  x  10  seconds,  four  additional  elements  experi¬ 
enced  increased  interlaminar  damage.  These  elements  are  shown  in 
figure  30.  The  last  remaining  interface  in  the  elements  was  the 
bottom  0°/45°  ply  junction.  At  the  end  of  the  analysis,  no  addi¬ 
tional  delamination  had  occurred. 

After  reviewing  the  results  of  the  second  stress  analysis,  it 
was  decided  that  a  third  should  be  performed.  The  decision  v/as  made 
in  an  effort  to  determine  the  damage  subsequent  to  the  point  at 
which  the  two  stress  solutions  had  terminated. 

Stress  Solution  with  Modified  Finite  Element  Model 

In  order  to  continue  the  solution  process  beyond  the  point  at 
which  an  element  has  no  remaining  plies,  two  approaches  were  con¬ 
sidered.  The  first  approach  involved  modifying  the  computer  code 
to  simply  ignore  the  element  after  the  damage  occurred.  The  time 
involved  in  adopting  this  method  was  not  available,  however,  and 
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a  simple  approach  was  needed.  The  approach  taken  involved  simply 
removing  the  four  central  elements  from  the  finite  element  model 
before  beginning  the  analysis.  It  was  felt  that  this  approach  would 
eliminate  the  most  critical  elements  and  allow  the  solution  to  pro¬ 
ceed  to  conclusion.  Since  these  four  elements  had  already  been 
shown  to  sustain  damage  throughout  the  laminate,  removing  them  would 
simulate  the  response  after  they  had  failed. 

The  finite  element  model  used  for  this  analysis  was  previously 
shown  in  figure  3  and  the  impact  mass  distribution  demonstrated  in 
figure  4 . 

The  solution  with  this  modified  finite  element  model  progressed 
through  30  time  steps,  as  had  the  original  stress  solution.  After 
30  steps,  the  solution  again  terminated  due  to  the  lack  of  material 
left  in  an  element.  The  progression  of  damage  growth  for  this  so¬ 
lution  is  shown  in  figures  31  through  37.  The  only  major  differ¬ 
ence  between  this  solution  and  the  previous  ones  relates  to  the 
location  of  the  element  which  caused  the  solution  to  terminate. 

In  the  previous  analyses,  the  element  which  failed  totally  had 
always  been  one  of  the  four  central  elements.  In  the  modified  model 
solution,  these  elements  were  removed.  The  element  which  caused  the 
execution  to  terminate  in  the  modified  model  solution  was  simply  an 
adjacent  element.  Hence,  the  effort  to  continue  the  analysis  be¬ 
yond  7.5  X  10”^  seconds  was  in  vain. 
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IV.  DISCUSSION  AND  RECOMMENDATIONS 


The  computerized  analysis  methodology  has  been  uti3„ized  to 
predict  the  response,  including  damage  initiation  and  propagation,-. 
for  a  specific  low  velocity,  low  mass  impact  event.  The  analysi.s 
has  been  -shown  to  be  an  effective  tool  in  predicting  the  response 
of  the  plate  structure.  The  damaged  nodes  and  locations  predicted 
were  consistent  and  correspond  reasonably  well  with  experimental 
work  as  described  to  MSC  by  NASA  Langley.  Several  significant 
features  of  the  analysis  predictions  and  methodology  warrant  fur¬ 
ther  discussions. 

0a^ch  of  the  stress  solutions  performed,  signij-icant  inter 
laminar  separations  occurred  at  the  first  stress  computation.  The 
time  at  which  the  calculations  were  made  vjas  1.25  x  10  seconds. 

The  elements  which  suffered  interlaminar  separation  at  this  time 
were  the  only  elements  affected  by  delamination.  These  two  features 
identify  two  significant  elements  of  the  stress  and  daraage  predic¬ 
tions  . 

The  fact  that  the  delamination  occurred  at  the  first  stress  solu¬ 
tion  demonstrates  that  interlaminar  failures  initiate  very  eaily  in 
the  impact  event.  The  delamination  is  a  function  of  the  transverse 
shear  forces,  which  are  highest  very  early  in  the  impact  event  be¬ 
cause  of  the  localized  nature  of  the  displacement  response  at  that 
time.  Since  the  significantly  non-zero  displacements  are  restricted 
to  a  small  region  at  the  center  of  the  plate,  severe  moment  gradients 
also  exist  at  the  plate  center.  The  high  moment  gradients  produce 
the  large  transverse  shear  forces.  As  the  impact  event  progresses, 
the  displacement  fields  spread,  and  the  moment  gradients  decrease 
and  the  propensity  for  interlaminar  damage  initiation  decreases. 

As  delamination  occurs,  shear  deformations  increase  due  to 
singular  shear  stresses  at  the  crack  tip.  The  singular  shear 
stresses  tend  to  promote  propagation  of  the  crack.  These  effects 
cannot  be  modeled,  however,  since  shear  deformation  is  not  in¬ 
cluded  in  the  current  analysis. 
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This  is  an  area  where  further  work  could  significantly  improve  the 
usefulness  of  the  CLIP  code.  The  propagation  of  interlaminar  dam¬ 
age  might  be  handled  in  one  of  two  ways. 

One  approach  would  involve  the  inclusion  of  shear  deformation 
in  the  displacement  solution.  With  the  inclusion  of  shear  defor¬ 
mation,  damaged  elements  could  be  given  a  reduced  stiffness  to 
promote  increased  shear  stresses  in  adjacent  elements  equal  to 
the  average  increase  produced  by  the  shear  stress  singularity. 

This  would  then  promote  delamination  propagation. 

The  other  approach  would  involve  the  use  of  an  analysis  to  be 
used  subsequent  to  the  impact  analysis.  This  method  would  use  the 
delamination  predicted  during  the  impact  analysis .  The  material 
surrounding  the  delamination  would  then  be  subjected  to  the  shear  and 
moment  distributions  predicted  in  the  impact  analysis  after  the 
instant  at  which  the  damage  occurred.  The  areas  where  the  delamina¬ 
tion  connected  with  intact  material  could  then  be  analyzed  utilizing 
fracture  mechanics  to  predict  the  delamination  growth.  Noting  that 
the  effects  of  shear  deformation  are  small  in  undamaged  materials, 
this  damage  growth  need  not  be  included  in  the  impact  analysis  if 
the  affected  area  is  small. 

Another  area  which  warrants  discussion  involves  the  extent  of 
the  ply  damage  predicted  in  the  impact  analyses  made.  In  all  of  the 
stress  solutions  made,  the  damage  eventually  propagated  completely 
through  the  thickness  of  the  central  finite  elements.  This  led  to 
an  immediate  problem  in  that  the  solution  terminated  at  this  point. 
Modifications  to  the  CLIP  code  to  allow  the  solution  to  proceed 
beyond  this  point  would  be  most  helpful .  The  solution  could  be  modi¬ 
fied  to  continue  until  the  damaged  region  reached  an  area  equivalent 
to  the  impacting  mass.  At  this  point  one  would  assume  that  the  mass 
had  penetrated  the  plate. 

Because  the  solutions  terminated,  the  full  extent  of  the  damage 
which  would  have  been  induced  in  the  plate  was  never  determined.  It 
was  apparent,  however,  that  the  predictions  of  damage  exceeded  the 
damage  measured  in  the  experimental  work  performed  at  NASA  Langley. 
This  was  primarily  related  to  the  damage  through  the  thickness 
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predicted  at  the  plate  center  and  this  would  tend  to  indicate  that 
the  energy  impacted  to  the  system  in  the  analysis  was  too  great. 

The  excess  energy  in  the  analysis  is  probably  a  function  of  two 
effects  not  modeled.  The  first  of  these,  and  probably  the  most 
iirportant,  is  related  to  energy  which  should  be  lost  due  to  local 
surface  crushing  at  the  point  of  impact.  It  is  known  from  Hertzian 
contact  analyses  that  the  local  contact  stresses  can  be  quite  large 
and  could  cause  significant  permanent  local  deformations.  These 
deformations  would  absorb  a  significant  portion  of  the  impact  energy. 
Modifications  to  account  for  this  effect  would  improve  the  analytical 
predictions.  These  modifications  could  be  effected  utilizing  the 
contact  forces  prediction  already  in  the  CLIP  code  coupled  with  a 
Hertzian  contact  analysis  and  a  local,  non-linear  material  model. 

The  other  contributor  to  the  energy  loss  is  material  damping. 

The  effects  of  damping  very  early  in  the  impact  event  are  probably 
small  but  at  later  times  they  must  surely  become  significant.  The 
CLIP  code  has  provisions  for  damping  but  this  feature  was  not  used 
due  to  a  lack  of  data. 

Another  area  of  discussion  relates  to  the  effects  of  utilizing 
a  distributed  impact  mass  in  the  analysis.  As  was  demonstrated, 
distributing  the  impact  mass  produced  a  slightly  increased  dis¬ 
placement  field  when  the  opposite  should  have  been  true.  The  effect 
was  seen  to  be  short-lived  and  may  be  insignificant.  It  is  discon 
certing,  however,  and  modifications  could  be  made  to  correct  the 
situation.  These  modifications  would  involve  applying  the  mass  and 
velocity  initially  at  one  node  and  spreading  the  mass  as  a  function 
of  the  contact  force  as  the  solution  progressed.  This  would  more 
closely  simulate  the  actual  impact  event. 

A  final  area  of  consideration  relates  to  the  effects  of  a 
static  prestress  on  the  impact  analysis .  The  procedure  utilized  in 
the  current  analysis  is  simply  a  superposition  of  the  static  and 
dynamic  displacement  fields.  This  is  applicable  if  in— plane  static 
loads  are  tensile  or  if  shear  or  compressive  loadings  are  small 
enough  that  buckling  is  not  a  consideration.  If  buckling  is  a  real 
possibility  then  the  inclusion  of  a  stability  analysis  would  be 
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required. 

A  stability  analysis  in  the  finite  element  procedure  would 
require  considerable  modifications  to  the  code  and  substantially 
decrease  the  efficiency  of  the  analysis.  A  stability  analysis 
requires  the  formulation  of  a  modified  stiffness  matrix  and  an 
eigenvalue  extraction.  The  modified  stiffness  matrix  is  based  on 
the  strain  field  present  when  the  buckling  analysis  is  performed. 
The  eigenvalue  problem  is  not  compatible  with  the  direct  integra¬ 
tion  of  the  equations  of  motion  as  currently  used  for  the  impact 
displacement  response.  Thus,  a  new  analysis  procedure  would  be 
needed.  Additionally,  the  eigenvalue  problem  would  have  to  be 
solved  at  each  time  step  since  the  buckling  loads  are  a  function  of 
the  strain  field  present.  This  would  effectively  double  the  solu¬ 
tion  time.  Thus,  while  the  inclusion  of  a  stability  analysis  would 
be  possible,  it  would  not  be  practical. 

In  svimmation,  while  the  analysis  is  functioning  well  and  pro¬ 
vides  significant  insight  into  the  phenomenon  of  low  velocity  im¬ 
pact,  further  work  would  be  most  helpful. 

The  areas  of  these  future  efforts  should  include: 

1)  Inclusion  of  a  disbond  growth  model  either 
through  shear  deformation  or  fracture  mechanics; 

2)  Modification  to  allow  continued  program  execution 
after  local  complete  element  failure; 

3)  Inclusion  of  energy  loss  mechanisms  due  to  non¬ 
linear  contact  effects;  and 

4)  More  realistic  modeling  of  the  impact  mass 
distribution . 

The  addition  of  these  capabilities  to  the  analysis  procedure 
would  provide  better  modeling  of  the  impact  event.  This  would  al¬ 
low  for  increased  confidence  in  the  analysis  and  facilitate  its  use 
in  evaluating  the  relative  merits  of  various  laminates  when  sub¬ 
jected  to  low  velocity  impact. 
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As  optimum  laminates  are  determined  from  impact  related  cri¬ 
teria,  further  evaluations  could  then  be  performed  with  regard  to 
residual  static  strength,  subsequent  impact  events,  and  fatigue 
life.  The  static  and  fatigue  evaluations  could  be  made  using  a 
modification  of  a  code  developed  for  fatigue  of  notched  composite 
laminates  (reference  18)  while  subsequent  impact  events  might  be 
evaluated  using  a  modification  of  the  code  developed  here. 
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Table  1.  Failure  Criteria 


Tensile  Fiber  Mode 


1 


Compressive  Fiber  Mode 


Tensile  In-Plane  Matrix  Mode  (<^22  ^  0*0) 


2  „2 
°22  ^  °12 
(a})  2  <t^>2 


1 


Compressive  In-Plane  Matrix  Mode 


(022  <  0.0) 


1]  + 


22 


12 


(2t^) 


1 
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Table  1  (cont'd.).  Failure  Criteria 


Tensile  Interlaminar  Mode  (<^33  ^  0.0) 


2  2  2 

®33  ,  ^23  ^13  _  , 

,  +.2'^,  ^2  ,  ^2 

(oj)  (t^)  (t^) 


Compressive  Interlaminar  Mode 


(0^3  <  0.0) 


1]  + 
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(2t^) 


(T^) 


^'^23‘*‘°13^ 


1 


where : 

=  axial  tensile  strength 
a~  =  axial  compressive  strength 
=  transverse  tensile  strength 
=  transverse  compressive  strength 
=  axial  shear  strength 
Ty  =  transverse  shear  strength 
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Table  2.  T300/5208  Properties 


”” 

153  GPa 

^T 

= 

10.9  GPa 

= 

5.6  GPa 

'’at 

= 

0.30 

P 

= 

1.55  g/cc 

0.70 

+ 

^A 

689.5  MPa 

758.5  MPa 

+ 

= 

27.6  MPa 

= 

96.5  MPa 

= 

62.1  MPa 

62.1  MPa 
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Table  3 .  Impact  Parameters 


IMPACT  MASS  =  16.45  g 

IMPACT  VELOCITY  =9.4  m/s 


Table  4. 


Comparison  Groupings  of  Predicted  Deunage 


Figures  14  - 

22 

Predicted  damage  in  solution  with  mass  lumped 
at  one  node 

Figures  23  - 

30 

Predicted  damage  in  solution  with  mass  dis¬ 
tributed  over  five  nodes 

Figures  31  - 

37 

Predicted  damage  in  solution  with  modified 
finite  element  model  (four  central  elements 
removed) 

ONE  NODE  IMPACT 

Figures 

Comparison 

14,  16  -  20 

Figures  show  the  extent  of  elements  with 
predicted  ply  damage  as  a  function  of  time 

14,  15 

A  comparison  of  the  prediction  of  elements 
with  ply  damage  vs.  elements  with  interlami¬ 
nar  damage  at  T  =  1.25  x  10”^  sec 

16,  17 

A  comparison  of  the  number  of  elements  with 
predicted  top  ply  damage  vs.  the  number  of 
elements  with  ply  damage  anywhere  through  the 
laminate  thickness  at  T  =2.5  x  10“^  sec 

20  -  22 

A  description  of  the  extent  of  ply  damage 
through  the  thickness  of  the  laminate 

FIVE  NODE  IMPACT 

Figures 

Comparison 

23,  25  “  27 

Figures  show  the  extent  of  elements  with 
predicted  ply  damage  as  a  function  of  time 

23,  24  A  comparison  of  the  prediction  of  elements 

with  ply  damage  vs.  elements  with  interlami¬ 
nar  damage  at  T  =  1.25  x  10“^  sec 

Figures  depict  the  extent  of  damage  in  the 
laminate  including  ply  damage,  fiber  damage 
and  interlaminar  damage 


27  -  30 


Table  4  (continued) .  Comparison  Groupings  of 
Predicted  Damage 


MODIFIED  FINITE  ELEMENT  MODEI. 

Figures  Comparison 

31,  33,  34  -  37  Figures  show  the  extent  of  elements  with  pre¬ 

dicted  ply  damage  as  a  function  of  time 

36,  37  A  comparison  of  elements  with  predicted  ply 

damage  vs.  elements  with  interlaminar  damage 
at  T  =  1.25  X  10~^  sec 


COMPARISONS  BETWEEN  ANALYSES 


Fic 

jures 

Comparison 

15, 

24, 

32 

A  comparison  of  the  number  of  elements  which 
have  sustained  interlaminar  damage  in  each 
of  the  three  solutions 

19, 

27, 

35 

A  comparison  of  the  number  of  elements  which 
have  sustained  ply  damage  at  T  =  5.0  x  10“^ 
sec  in  each  of  the  three  solutions 

20, 

27, 

37 

A  comparison  of  the  number  of  elements  with 
predicted  ply  damage  at  the  end  of  each  solu¬ 
tion 

20, 

28 

A  comparison  of  the  number  of  elements  which 
have  sustained  fiber  damage  in  the  solution 
with  one  impacted  node  vs.  the  solution  with 

five  impacted  nodes 
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Figure  3.  Modified  Finite  Element  Model 
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Figure  5.  Displacement  Response  Through  the  Plate  center 
Along  the  Shorter  Axis,  Displacement  Solution, 
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Figure  6.  Displacement  Response  Through  the  Plate  Center 
Along  the  Shorter  Axis  ,  Displacement  Solution 


Figure  12. 


Transverse  Shear  Resultant,  Q,  Through  the  Plate 
Center  Along  the  Shorter  Axis^  Stress  Solution, 
Mass  Lximped  at  One  Node  . 
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Figure  13.  Transverse  Shear  Resultant,  Q  ,  Through 
the  Plate  Center  Along  the  Longer  Axis, 
Stress  Solution,  Mass  Lumped  at  One  Node 


Elements  With  Ply  Damage,  T  =  1.25  x  10 
sec..  Stress  Solution,  Mass  Lumped  at  One  Node 
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Figure  16.  Elements  with  Ply  Damage,  T  =  2.5  x  10~^ 

sec..  Stress  Solution,  Mass  Lumped  at  One  Node 
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Figure  18-  Elements  with  Ply  Damage,  T  =  3.75  x  10 

sec..  Stress  Solution,  Mass  Lumped  at  One  Node 
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Figure  20.  Elements  with  Ply  Damage,  T  =7.5  x  10~^  sec 
Stress  Solution,  Mass  Lumped  at  One  Node 
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Figure  23.  Elements  with  Ply  Damage,  T  =  1.25  x  10”^ 

sec..  Stress  Solution,  Mass  Distributed  at  5 
Nodes 
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Elements  with  Ply  Damage,  T  =  3.75  x  10 
sec.,  Stress  Solution,  Modified  Model 
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Figure  35.  Elements  with  Ply  Damage,  T  =  5.0  x  lo"*^ 
sec..  Stress  Solution,  Modified  Model 


APPENDIX  A 


STRESS,  MOMENT,  AND  TRANSVERSE  SHEAR  RESULTANTS 

In  order  to  perform  the  composite  stress  analysis  required  for 
failure  prediction,  it  is  necessary  to  determine  stress,  moment, 
and  transverse  shear  resultants.  The  thin  shell  finite  element 
taken  from  SAP  IV  has  the  capability  for  stress  and  moment  result¬ 
ants  but  lacks  the  ability  to  compute  transverse  shear  resultants . 

It  was  necessary,  therefore,  to  add  the  capability  for  this  trans¬ 
verse  shear  prediction. 

Since  the  primary  focus  of  the  current  contract  involved  ma¬ 
terial  response  and  not  the  development  of  general  purpose  finite 
element  routines,  the  simplest  approach  available  for  computing 
transverse  shear  resultants  was  adopted.  In  the  course  of  this 
effort,  it  was  determined  that  since  the  restrictions  inherent  in 
the  simple  transverse  shear  computations  effectively  remove  the  ne¬ 
cessity  for  the  con^lex  stress  and  moment  resultant  computation, 
the  SAP  IV  routines  which  perform  these  analyses  were  not  required. 
Therefore,  the  entire  stress  recovery  procedure  is  based  on  the 
simplifying  element  orientations  required  for  transverse  shear 
as  described  below. 

The  primary  restriction  invoked  by  the  simple  transverse  shear 
resultant  computation  involves  the  elemental  geometry.  No  attempt 
was  made  to  determine  transformations  of  forces  from  the  global 
coordinate  system  to  an  arbitrary  elemental  coordinate  system. 
Therefore,  it  is  required  that  the  local  elemental  coordinates  coin¬ 
cide  with  the  global  coordinates.  Additionally,  no  provisions  were 
made  for  elemental  geometric  irregularities.  This  requires  that  all 
elements  be  rectangular.  The  second  restriction  applies  only  to  the 
stress  recovery  procedures.  The  displacement  response  is  computed 
correctly  for  elements  which  are  not  rectangular  if  the  element  co¬ 
ordinate  system  is  defined  coincident  with  the  global  coordinate 
system.  This  restriction  is  required  so  that  the  laminate  prop¬ 
erties  are  utilized  correctly. 

In  figures  A-1  and  A-2,  the  elemental  nodal  forces  and  moments 
are  depicted.  The  forces  and  moments  are  shown  in  a  global  positive 
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scene.  The  element  is  required  to  be  rectangular  and  oriented  with 
the  global  coordinate  system  as  mentioned  previously  and,  thus,  the 
relations  for  stress  and  moment  resultants  can  be  written  directly 
as  given  in  tables  A-1  and  A-2.  It  can  easily  be  seen, that  the 
relations  in  tables  A-1  eind  A-2  are  invalid  if  either  of  the  two 
restrictions  mentioned  are  violated. 

The  computation  of  transverse  shear  stress  resultants!  is  slight¬ 
ly  more  complicated.  In  figvire  A-3  the  moments  and  shear  resultants 
required  for  equilibrium  are  depicted.  In  the  analysis  used,  the 
transverse  shear  forces  are  computed  from  moment  equilibrium  since 
shear  deformation  is  not  included  in  the  displacement  solution. 

The  transverse  shear  resultants  can  be  determined  quite  easily 
from  equilibrium  considerations  as: 


9M  3M 

Q  =  -  _JQ1  _ 

y  9x  3y 


(A^l) 


If  it  is  assumed  that  all  quantities  vary  linearly  within  the 
finite  element,  equations  [A.l]  can  be  rewritten  as: 


Each  of  the  quotients  in  equations  [A. 2]  can  then  be  written  in 
terms  of  the  nodal  moments  of  figure  A-2.  When  this  is  done,  a  dif¬ 
ficulty  is  immediately  encountered. 
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The  various  terms  of  equation  [A. 2]  are  found  to  be: 
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(A. 3) 


It  is  immediately  apparent  that  in  terns  of  the  nodal  moments, 
the  two  terms  required  for  are  identical,  as  are  the  two  terms 
in  the  expression.  This  result  is  demonstrating  that,  in  terms 
of  the  nodal  moments,  it  is  not  possible  to  separate  the  contribu¬ 
tions  of  the  bending  and  twisting  moments.  This  is  easily  confirmed 
by  considering  a  one-dimensional  case.  The  shear  force  computed  is 
twice  the  proper  value  if  the  full  expressions  of  equations  [A. 2] 
are  utilized. 

The  resolution  of  this  problem  involves  simply  dropping  the 
twisting  components  from  equation  [A. 2]  and  using  the  expressions: 


AM 

- - Z 

Ay 


(A. 4) 
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It  must  be  remembered  that  the  contributions  of  the  twisting 
moments  have  not  been  neglected  in  equations  [A.  4] .  These  contribu' 
tions  are  included  implicitly  in  equations  [A. 4]  through  the  confu¬ 
tations  of  the  two  quotients  involved. 


Table  A-1.  Stress  Resultants 
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Table  A-2 


Moment  Resultants 
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APPENDIX  B 


INTERLAMINAR  SHEAR  STRESSES 

The  analysis  used  in  the  CLIP  code  for  computing  interlaminar 
shear  stresses  is  based  upon  bending  equilibrium.  In  Appendix  A, 
the  derivation  of  the  transverse  shear  stress  resultants  is  given. 
These  resultants  are  then  used  to  compute  transverse  shear  stresses 
at  the  ply  interfaces  within  the  laminate.  The  method  used  for  com¬ 
puting  the  interlaminar  shear  stresses  is  based  on  an  extension  of 
classical  methods  used  for  bending  shear  stresses  in  homogeneous 
beams. 

In  the  analysis,  it  is  necessary  to  compute  interlaminar  shear 
for  two  different  cases.  Interlaminar  stresses  must  be  computed  in 
elements  which  have  sustained  damage  as  well  as  those  which  have 
not.  The  methods  used  in  both  cases  are  described  here. 

INTERLAMINAR  SHEAR  IN  INTACT  ELEMENTS 

In  order  to  compute  interlaminar  shear  stresses  in  undamaged 
elements,  the  transverse  shear  stress  resultants  are  converted  to 
moments  by  multiplying  by  the  appropriate  shell  element  dimensions 
(see  eqn.  A. 4).  These  two  moments  then  represent  the  change  in 
moment  across  the  element  in  both  the  X  and  Y  coordinates .  This  is 
demonstrated  for  a  one-dimensional  case  in  figure  B-1.  It  is  readily 
apparent  in  figure  B-1  that  the  shear  stresses  represent  the  balanc¬ 
ing  force  required  for  moment  equilibrium. 

The  moment  differentials  derived  from  the  transverse  shear 
resultants  are  then  applied  to  the  laminate  model.  This  produces 
three  in-plane  stresses  within  each  ply.  These  ply  stresses  are 
then  converted  to  forces  and  summed  through  the  thickness  in  both 
the  X  and  Y  coordinate  directions .  In  figure  B-2 ,  the  three  in¬ 
plane  stresses  and  resulting  interlaminar  shear  stresses  are  shown 
for  an  outer  ply  of  the  laminate  and  an  applied  moment  differential. 
Equilibrium  in  the  X-direction  yields  forces 
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Ay  Az  +  Az  +  Ax  Ay  =  0 


(B.l) 


and  in  the  Y-direction 

Oy  Ax  Az  +  0-^  Ay  Az  +  Ax  Ay  =  0. 

Thus,  the  interlaminar  shear  stresses  are  simply 

-  (o  Ay  Az  +  a  Ax  Az) 

X  xy 

Ax  Ay 


(B.2) 


{B.3) 


To  compute  interlaminar  shear  stresses  on  interior  ply  inter¬ 
faces,  it  is  only  required  that  the  X  and  Y  forces  represented  by 

the  0  and  o  terms  respectively  in  equation  B.l  be  summed  through 
■xy  yz 

the  laminate  to  the  interface  in  question. 


INTERLAMINAR  SHEAR  IN  DAMAGED  ELEMENTS 

In  order  to  compute  interlaminar  shear  stresses  in  an  ele¬ 
ment  which  has  sustained  damage  prior  to  the  stress  calculation, 
it  is  first  necessary  to  determine  how  the  laminate  responds  to 
loading  when  damage  is  present.  In  figure  B-3  the  bending 
stresses  through  a  localized  delamination  are  shown.  Since  the 
delamination  is  small,  the  material  away  from  the  delamination 
effectively  forces  the  curvature  and  strain  field  to  remain  un¬ 
changed.  The  material  above  and  below  the  delamination  bend  as 
two  independent  laminates  while  the  constraint  of  the  adjacent, 
undamaged  material  adds  opposing  membrane  forces.  Thus,  the  net 
effect  produces  no  change  in  the  bending  stress  field.  This  is 
the  basis  for  not  reformulating  the  elemental  stiffness  for 
interlaminar  delaminations. 
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This  model  cannot  account  for  a  moment  gradient,  however.  In 
order  to  evaluate  the  effects  of  a  moment  gradient,  it  is  required 
that  the  two  sublaminates  behave  entirely  independently .  The  op¬ 
posing  membrcine  forces  must  not  be  present  in  the  shear  analysis. 

If  a  moment  gradient  is  applied  to  this  type  of  model,  then  the  mem¬ 
brane  forces  must  also  exhibit  a  gradient.  This  is  not  possible 
since  it  would  require  a  shear  force  transfer  across  the  delamination. 
This  is  a  drawback  related  to  the  lack  of  shear  deformation  in  the 
analysis. 

To  model  the  moment  gradient,  it  is  necessary  to  compute  the 
bending  stiffnesses  of  the  sublaminates  independently.  The  bending 
stiffnesses  are  then  added  to  give  the  stiffnesses  of  the  assemblage 
of  partial  laminates.  This,  however,  leads  to  bending/extensional 
coupling  at  both  the  laminate  and  sublaminate  levels.  It  is  required 
therefore,  that  the  [B]  matrices  be  accounted  for.  This  can 
easily  be  accomplished  by  noting  that  the  membrane  forces  must  be 
zero  for  the  applied  transverse  shear  forces  (moment  gradients). 

The  force/moment  -  strain/curvature  relations  for  a  general 
laminate  are 


{n}  =  [A]{e°}  +  [B]{ic} 
{m}  =  [B]{e°}  +  [D]{ic} 


(B.4) 


where  the  matrices  have  the  usual  connotations.  Remembering  that 


{n}  =  0 


(B.5) 


the  mid-plane  strains  due  to  a  bending  load  are 

{e°}  =  (B.6) 

These  strains  actually  represent  a  neutral  axis  shift  from  the 
mid-surface  of  the  unsymmetric  laminate  to  the  proper  neutral  surface 
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Using  the  strains  of  equation  (B.6), 
stiffness  matrix  can  be  determined. 

an  effective  bending 

{m}  =  [B]{e°}  +  [D]{tc} 

(B.7) 

=  -IB]  [A]"^  [B]{k}  +  [D]{tc} 

(B.8) 

=  ( [D]  -  [B] IA]“^  IB])  {<} 

(B.9) 

=  ID*]{k} 

(B.IO) 

The  bending  stiffness  matrix  in  (B.IO)  relates  moments  and  curva¬ 
tures  about  the  neutral  surfaces  rather  than  the  mid-surface. 

The  bending  stiffness  matrix  for  the  total  laminate  is  then 
given  as 

'  .N  ,  "  ^ 

ID*]  =  Z  [D*]  ,  N  -  number  of  sublaminates.  (B.ll) 

■V',,'  ,  :  1=1  ■ 

with  the  bending  stiffness  of  the  assemblage  of  sublaminates  given 
in  egn.  (B.ll) ,  a  curvature  vector  can  easily  be  found  for  the  ap¬ 
plied  moment  differentials.  Using  the  curvature  for  the  whole  lami¬ 
nate  and  eqn.  (B.6) ,  the  ply  stresses  for.  each  sublaminate  can  be 
determined  and  summation  of  forces  depicted  in  figure  B-2  and  equa¬ 
tions  B.l  through  B.3  carried  out.  Thus,  the  interlaminar  stresses 
in  damaged  elements  are  predicted. 

When  an  element  contains  both  ply  damage  and  interlaminar  damage 
the  same  method  is  used.  In  figure  B-4  the  sublaminates  for  a  lami¬ 
nate  which  has  sustained  general  damage  are  shown.  For  ply  damage, 
the  entire  ply  is  removed  from  interlaminar  shear  calculations. 
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Figure  B-4.  Typical  Failed  Laminate 
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APPENDIX  C 


COMPOSITE  LAMINATE  IMPACT  PROGRAM  (CLIP) 

The  computer  program  developed  for  predicting  damage  initiation 
and  growth  during  low  velocity  impact  consists  of  a  transient  dynam¬ 
ics  finite  element  code  coupled  with  composite  stress  cind  failure 
analysis  procedures.  In  conjunction  with  the  failure  analysis  rou¬ 
tines  is  the  capability  to  incorporate  predicted  damage  into  the 
trainsient  analysis.  Thus,  the  damage  predicted  during  any  time  step 
is  incorporated  into  the  dynamic  solution  at  future  time  steps. 

This  coupling  of  damage  and  dynamic  response  is  the  heart  of  the  com¬ 
puterized  procedure. 

The  transient  dynamic  finite  element  procedures  used  were  taken 
from  the  SAP  IV  code  (ref.  14).  The  portions  of  the  SAP  IV  code 
used  include  the  thin  shell  element  routines  and  the  time  integra¬ 
tion  routines.  The  composite  stress  and  failure  analysis  routines 
were  developed  specifically  for  use  in  the  CLIP  code. 

In  figure  G-1,  a’  simplified  flow  chart  of  the  analysis  performed 
by  the  CLIP  code  is  shown.  The  branches  and  loops  shown  represent 
the  procedures  required  for  incorporating  the  effects  of  predicted 
damage  into  the  analysis  and  removing  the  impact  mass  after  the 
contact  force  becomes  tensile.  The  various  computer  program  seg¬ 
ments  which  perform  the  analysis  are  listed  in  table  C-1  along  with 
a  brief  description  of  their  functions. 

PROGRAM  OPERATION 

The  CLIP  code  has  been  written  with  several  different  solution 
options.  The  different  options  offer  the  user  considerable  flexi¬ 
bility  in  using  the  code.  The  options  available  include  data  check 
mode,  static  analysis  alone,  static  analysis  to  generate  a  pre¬ 
stress  condition  for  a  dynamic  analysis ,  dynamic  analysis  without 
stress  calculation  and  dynamic  analysis  with  stress  calculation. 

The  data  check  mode  of  operation  is  useful  in  insuring  that 
all  input  data  parameters  are  correct.  This  option  is  essential 
in  any  large  program. 
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The  static  solution  mode  has,  basically,  two  uses.  First,  a 
static  test  case  is  an  easy  method  to  prove  the  validity  of  a  finite 
element  model.  It  is  extremely  difficult  to  verify  the  results  of 
a  dynamic  analysis  directly.  The  second  use  of  the  static  analysis 
is  to  generate  a  pre-stress  condition  at  the  beginning  of  the  dynamic 
analysis.  This  allows  the  static  stresses  to  be  included  directly 
with  the  dynamically  induced  stresses. 

The  inclusion  of  the  static  solution  is  performed  by  direct 
superposition.  The  static  displacement  vector  is  saved  at  the  end 
of  each  time  step.  In  this  way  the  static  displacements  do  not  act 
as  an  initial  condition  in  the  dynamic  response.  It  must  be  remem¬ 
bered  that  this  is  a  direct  superposition  of  results.  There  is  no 
coupling  of  the  static  and  dynamic  responses. 

The  dynamic  solution  can  be  performed  with  or  without  stress 
calculations.  The  dynamic  analysis  made  without  stress  calculations 
is  considerably  less  costly  and  time  consuming  than  a  similar  analy¬ 
sis  with  a  complete  stress  analysis.  The  solution  cost  reductions 
are  a  result  of  two  factors.  First,  the  stress  analysis  procedures 
are  by-passed  in  the  analysis  and  secondly,  there  are  no  damage  cal¬ 
culations  which  would  require  reformulation  of  elemental  and  global 
stiffness  matrices  and  subsequent  decomposition.  Because  of  this, 
the  dynamic  solution  without  stress  analysis  can  be  used  to  determine 
the  length  of  time  required  to  characterize  the  impact  event  at  a 
moderate  cost.  This  mode  of  analysis  can  also  be  used  to  verify 
the  integration  time  step  selected. 

When  using  the  dynamic  solution  without  stress  analysis,  it 
should  be  remembered  that  the  response  of  the  plate  will  be  different 
when  stress  calculations  are  included.  The  amount  of  this  difference 
will  be  a  function  of  the  amount  of  damage  sustained  in  the  plate 
when  stress  calculations  are  included. 

When  damage  predicted  in  an  element  is  ply  damage  (either  matrix 
or  fiber) ,  the  elemental  stiffness  matrix  is  modified.  These  modi¬ 
fications  are  then  included  in  the  global  stiffness  matrix  for  in¬ 
corporation  in  subsequent  time  steps.  Thus,  the  plate  stiffness  is 
changed  forcing  a  plate  response  change. 
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The  dynamic  solution  mode  which  includes  stress  analysis  is  the 
primary  mode  of  operation  of  the  ctilP  code.  This  mode  is  used  for 
simulating  the  actual  response  of  a  composite  plate  subjected  to  low 
velocity  impact. 

The  use  of  the  analysis  in  this  mode  requires  some  thought  on 
the  part  of  the  user.  The  code  is  written  such  that  the  user  selects 
how  often  the  stress  calculations  are  performed.  If  the  stress  cal¬ 
culations  are  not  performed  often  enough,  it  is  possible  to  miss 
considerable  damage.  This  is  especially  true  at  the  beginning  of 
the  analysis  when  very  localized  deformations  are  the  predominant 
response  of  the  plate. 

In  an  effort  to  achieve  computational  economy,  the  CLIP  code 
performs  stress  analysis  only  for  elements  which  the  user  selects. 
This  option  saves  considerable  cbmputei'  time  since  many  of  the 
finite  elements  in  the  model  will  typically  never  sustain  daunage. 

If,  however,  an  element,  which  should  fail  under  the  dynamic  load¬ 
ing,  is  not  selected  for  stress  analysis,  this  damage  is  lost.  It 
is  never  computed  and,  hence,  is  not  included  in  the  dynamic  response 
Thus,  it  is  important  to  select  the  proper  group  of  elements  for 
stress  calculations.  Insight  into  the  elements  which  should  be  in¬ 
cluded  can  be  obtained  by  excimining  the  plate  response  without  Stress 
calculations. 

PROGRAM  OUTPUT 

The  output  of  the  CLIP  code  can  be  divided  into  two  distinct 
entities.  First,  the  program  echos  the  input  data  such  that  a 
record  of  this  data  is  generated,  and  secondly,  the  results  of  the 
analysis  are  printed.  In  each  of  these  two  areas,  the  user  has 
considerable  flexibility  in  determining  the  extent  of  the  output 
generated. 

In  terms  of  the  input  data  echo,  the  user  can  suppress  printing 
of  the  nodal  points  or  the  elemental  data  or  both .  Other  input  data 
including  the  impact  parameters,  lamina  materials  and  laminate  con- 
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figuration  as  well  as  data  generated  by  the  analysis  code  with  re¬ 
spect  to  problem  size  and  laminate  elastic  constants  are  always 
printed. 

In  terms  of  the  solution  output,  the  user  selects  the  frequency 
of  output  for  the  various  data  generated  as  well  as  the  nodal  points 
(displacements,  velocities,  accelerations)  and  elements  (element 
forces,  laminate  coordinate  stresses,  layer  coordinate  stresses, 
failure  calculations,  failure  locations)  where  data  are  to  be  printed. 

For  the  nodal  output  data,  the  user  can  select  different  print¬ 
out  frequencies  for  the  displacements,  velocities  and  accelerations 
as  long  as  the  velocity  and  acceleration  print  frequencies  are  even 
multiples  of  the  displacement  frequency.  Thus,  the  displacements 
may  be  printed  every  five  time  steps  while  the  velocities  print 
every  fifty  steps  and  the  accelerations  every  20  steps.  The  print 
frequencies  are  left  to  the  discretion  of  the  user. 

The  elemental  data  are  printed  in  a  similar  fashion  where  each 
printout  frequency  must  be  an  even  multiple  of  the  stress  calculation 
frequency.  As  an  example,  the  user  must  not  select  a  stress  calcula¬ 
tion  frequency  of  every  ten  steps  while  requesting  a  layer  coordinate 
stress  printing  every  three  time  steps.  The  result  of  this  particu¬ 
lar  arrangement  would  produce  layer  stress  printing  only  every  thirty 
time  steps.  There  is  no  requirement  that  the  individual  data  print 
frequencies  be  multiples  of  each  other,  however. 

The  elemental  data  which  can  be  printed  consist  of  elemental 
forces,  ply  stresses  in  two  coordinate  systems,  calculations  of  the 
failure  analyses,  and  locations  of  damage.  The  elemental  forces 
consist  of  stress,  moment  and  transverse  shear  resultants  and  top 
surface  stress  corresponding  to  the  impactor/plate  contact  force 
distributed  over  all  elements  adjoining  the  impacted  nodes. 

Layer  stress  can  be  printed  in  the  laminate  coordinate  system 
(global  X-Y)  or  the  local  ply  coordinate  systems,  or  both.  Stresses 
are  printed  for  each  ply  in  the  laminate  under  study  as  well  as  each 
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ply  interface  including  the  top  and  bottom  free  surfaces.  Thus, 
there  is  always  one  more  interface  than  there  are  plies. 

The  calculations  of  the  failure  criteria,  if  selected  for 
printing,  produce  data  for  each  ply  and  interface.  The  data 
printed  are  simply  the  sum  of  the  terms  in  the  various  failure 
criteria  modes  (see  table  1) .  These  data  can  be  used  to  iden¬ 
tify  locations  and  modes  of  danage  within  the  elemental  laminate 
models. 

If  the  details  of  these  calculations  are  not  required,  the 
user  can  select  that  failure  location  printing  be  activated. 

This  output  details  the  ply  or  interface  which  has  failed  within 
a  damaged  element.  For  the  case  of  ply  damage,  information  is 
also  printed  indicating  whether  matrix  or  fiber  failure  has 
occurred.  It  is  possible  to  specify  both  failure  calculation 
printing  and  failure  location  printing.  To  do  so,  hpwevef,  is 
somewhat  redundant. 
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CLIP  USERS  GUIDE 


Composite  Laminate  Impact  Program 
I.  PROGRAM  CONTROL  DATA 


Card  1 

Title  Card  20A4 

Columns 

Contents 

1-80 

RED (20) 

Program  Title  Card 

Card  2 

Control  Data  415 

Columns 

Contents 

1-5 

NUMNP 

Number  of  nodes 

6-10 

MODEX 

Execution  code,  =0  execute 

11-15 

NSTR 

Stress  calculation  code 

-0,  no  stress  calculation 
=N,  stress  Calculation  every  N 
time  steps 

16-20 

ISTAT 

Static  analysis  code 
=0,  no  static  pre-stress 

Card  3 

Dynamic  Analysis  Control  Data  I5,3F10.0 

Columns 

Contents 

1-5 

NT 

Maximum  number  of  time  steps 

6-15 

DT 

Time  step 

16-25 

ALFA 

Mass  proportional  damping 

coefficient 

26-35 

BETA 

Stiffness  proportional  damping 

coefficient 
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Card  4 


Print  Control  Data 


1015 


Columns 

Contents  . 

1-5 

KEY(l) 

Node  print  code> 

^0  print  node  data  :  . 

6-10 

KEY (2) 

Element  print  code,  /O  print 

element  data 

11-15 

KEY (3) 

Displacement  print  code 

16-20 

KEY(4) 

Velocity  print  code 

21-25 

KEY(5) 

Acceleration  print  code 

26-30 

KEY(6) 

Element  force  print  code 

30-35 

KEY (7) 

Laminate  coordinate  stress  . 

pr  int  code 

36-40 

KEY (8) 

Lamina  coordinate  stress 

print  code 

41-45 

KEY (9) 

Failure  calculation  print 

code 

46-50 

KEY (10) 

Failure  location  print  code 

NOTES; 

Card  4 

KEY (3) -KEY (10) 

=0,  suppress  printing 

=N,  print  every  N  time  steps 

KEY (6) -KEY (10) 

Codes 

not  used  if  NSTR,  Card  2,  =0 

Card  5 

Print 

Control  Data  2015 

Column 

Contents 

1-5 

NDISP 

Number  of  nodal  print  groups. 

6-10 

NSTRP 

Number  of  element  stress 

calculation  groups  £100 
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NOTES ;  Card  5 

If  NDISP  =0,  print  data  for  all  nodes, 
do  not  include  Card  6 

If  NSTRP  =0,  compute  stresses  for  all  elements, 
do  not  include  card  7 


Card  6 

Print  Control  Data 

2015 

Column 

Contents 

1-5 

KEYS 

(1,1) 

First  node,  print 

group  I 

6-10 

KEYS 

(1,2) 

Last  node,  print 

group  I 

11-15 

KEYS 

(J,l) 

First  node,  print 

group  J 

16-20 

KEYS 

(J,2) 

Last  node,  print 

group  J 

Continue  through  NDISP  groups,  more  than  1  card 
if  necessary  , 


Card  7 

Print  Control  Data 

2015 

Column 

Contents 

1-5 

KEYS 

(1,1) 

First  element 

stress  calculation. 

Group  I 

6-10 

KEYS 

(1,2) 

Last  element  ! 

stress  calculation. 

Group  I 

11-15 

KEYS 

(J,l) 

First  element 

stress  calculation. 

Group  J 

16-20 

KEYS 

(J,2) 

Last  element. 

print  group  J 

Continue  through  NSTRP  groups,  more  than  1  card 
if  necessary 

NOTES ;  Cards  5,7 

Stress  calculations  are  performed  only  for  element 
identified  in  these  groups. 
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II.  IMPACT  CONTROL  DATA 


Card  1 


Node  Data 


15 


Columns 

1-5  NDMIMP 


Contents 

Number  of  impacted  nodes 
0  <  NDMIMP  <  20 


Card  2 


Node  Data 


2015 


Columns 

1-5  IMN(I) 

6-10  IMN(I+1) 


Contents 
Impacted  node 
Impacted  node 


Continue  through  NDMIMP  nodes. 


Card  3 

Impact 

Conditions 

2F10.0 

Columns 

Contents 

1-10 

TOTMAS 

Total  impactor  mass 

11-20 

TERMV 

Impact  velocity 

Card  4 

Impact 

Mass  Distribution 

8F10.0 

Columns 

Contents 

1-10 

XMFRAC(I) 

Impact  mass  fraction. 

node  I 

11-20 

XMPRAC(J) 

Impact  mass  fraction. 

node  J 

Continue  through  NDMIMP  nodes. 
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III.  LAMINATE  DATA 


Card 

1 

Laminate  Geometry  Control  Data  215 

Columns 

Contents 

1-5 

NMAT 

Number  of  different 

materials  £  2 

6-10 

NPLY 

Number  of  plies  £  25 

Card 

2 

Lamina  Elastic  Constants  5F10.0 

Columns 

I- 10  El (I) 

II- 20  E2(I) 

21-30  G(I) 

31-40  ANU(I) 

41-50  RHO{I) 

Continue  through  NMAT  itiaterials. 

Notes  ;  Card  2 

Major  Poissons  Ratio,  Vi2  =  ^ll'^2l/^22 
Card  3  Lamina  Material  Strengths  .  6F10.0 

Columns 

I- 10  STREN(1,I) 

II- 20  STREN(2,I) 


Contents 

Axial  tensile  strength,  material  I 
Axial  compressive  strength, 
material  I 


Contents 

Lamina  Axial  Modulus  En, 
material  I 

Lamina  Transverse  Modulus, 

E22f  material  I; 

Lamina  Axial  Shear  Modulus, 

Gi2r  material  I 

Lamina  Major  Poissons  Ratio 

■^12,  material  I 

Lamina  Mass  Density,  material  I 


104 


21-30 

STREN(3,I) 

Transverse 

tensile  strength. 

material  I 

31-40 

STREN(4,I) 

Transverse 

compressive  strength 

material  I 

41-50 

STREN(5,I) 

Axial  shear  strength,  material  : 

51-60 

STEEN (6, I) 

Transverse 

shear  strength ,  ^ 

material  1 

Continue 

through  NMAT  materials. 

Card  4 

Lamina  Geometries 

2f10.0,I5 

Contents 
Ply  thickness 
Ply  orientation,  degrees 
Ply  material  number 

Continue  through  NPLY  plies. 


Columns 

I- 10  THIK(I) 

II- 20  THET(I) 

21-25  MAT (I) 
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IV.  NODAL  DATA 


Card  1 

Columns 

I- 5 
6-10 

II- 15 

16-20 

21-25 

26-30 

31-35 

36-45 

46-55 

56-65 

NOTES ; 


Nodal  Input  Data 


7I5,3F10.0 


N 

ID(N,1) 

ID(N,2) 

ID(N,3) 

ID(N,4) 

ID(N,5) 

ID(N,6) 

X(N) 

Y(N) 

Z(N) 


Contents 

Nodal  point  number 
Boundary  condition  code, 
global  X  -  direction 
Boundary  condition  code, 
global  y  -  direction 
Boundary  condition  code, 
global  Z  -  direction 
Boundary  condition  code, 
global  X  -  rotation 
Boundary  condition  code 
global  y  -  rotation 
Boundary  condition  code 
global  Z-  rotation 
X  -  Coordinate 
y  -  Coordinate 
Z  -  Coordinate 


Continue  through  NUMNP  nodes 


Card  1, 

ID(I,J)=0,  Force  boundary  condition 

=1,  Zero  displacement  boundary  condition 
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V.  ELEMENT  DATA 


Card  1 

Columns 

1-5 

Card  2 

Columns 

I- 5 
6-10 

II- 15 
16-20 
21-25 
26-30 


31-40 


NOTES : 


Element  Input  Data  15 

Contents 

NUMEL  Number  of  elements 

Element  input  data  6I5,F10.0 


Contents 


MM 

Element 

number 

IY(1) 

Element 

node  I  - 

IY(2) 

Element 

node  J 

IY(3) 

Element 

node  K 

IY(4) 

Element 

node  L 

IY(5) 

Element 

stiffness  re-use 

code 

=0,  form  new  stiffness 
“1,  re-use  previous  stiffness 
PRESSU  Uniform  lateral  pressure 

load 

Continue  through  NUMEL  cards.  . 

Card  2 

IY(1)  -  IY{4)  , 

1.  Element  must  be  rectangular 

2.  Nodes  IY(1),  IY(2)  must  lie  along  global  X-Axis 

3.  IY(1)“IY{2)  Define  local  positive  X-Axis 

4.  Element  must  be  defined  counter-clockwise 

IY(5),  for  stiffness  matrix  re-use. 

Element  geometry  must  be  identical  to 
previous  element 
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PRESSU,  Used  in  static  analysis  only 

Ignored  if  ISTAT  =0,  Section  I,  Card  2 
VI  STATIC  NODAL  LOADS 

Card  1  Nodal  Load  Input  I5,F10.0 


Columns 
1-5  N 

6-15  TR(1) 

16-25  TR(2) 

26-35  TR(3) 

36-45  TR(4) 

46-55  TR(5) 


Contents 
Node  number 
X-Axis  load 
Y-Axis  load 
Z-Axis  load 
X-Axis  moment 
Y-Axis  moment 


NOTES ;  Card  2 

Used  in  static  analysis  only,  but  at  least  one  blank 
card  must  be  included. 


N  -  If  N*0  or  blank,  terminate  nodal  load  input 


Z-  Axis  moments  are  not  permitted  since  only  flat 

models  are  allowed  and  thus  no  Z-rotational  stiffness 
exists 


-  END  OF  INPUT  DATA  - 
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GENERAL  PROGRAM  NOTES: 


1.  Units  for  the  various  input  parameters  need  only  be  consistent 
with  the  exception  of  ply  orientation  angles,  which  must  be 
input  as  degrees. 

2.  When  inputting  nodal  boundary  condition  codes,  it  is  advis¬ 
able  to  use  zero-displacement  wherever  possible.  These 
serve  to  reduce  program  size  and  rimning  time.  Specifically, 
Z-axis  rotations  should  be  excluded.  X  and  Y  axis  displace¬ 
ments  should  be  excluded  unless  membrane  forces  are  included 
in  static  prestress  conditions. 

3.  Z  -  coordinate  location  must  be  identical  for  all  nodes. 

This  restriction,  along  with  the  restriction. that  elements  be 
rectangular  and  aligned  with  the  global  X-axis,  is  required 
for  proper  stress  calculation. 

4 .  Damping  coefficients  need  not  be  included  if  no  damping  is 
desired  in  the  dynamic  analysis . 

5.  The  laminate  input.  Section  III,  should  be  symmetric  as  the 
displacement  solution  ignores  bending/extensional  coupling. 


TABLE  C-1.  CLIP  ROUTINES 


CLIP 


TIMER 

IMPIN 

LAMIN 


INPUTJ 


ELT6 

TPLATE 

STRETR, 

QTSHEL, 

QOCOS , 

TDCOS, 

TRFPRD, 

SLST, 

SLCCT 

CALBAN 


Main  Program.  Supervises  the  data  input.  Sets 
internal  storage  parameters.  Supervises  element 
and  global  stiffness  formulation. 

Computes  and  prints  elapsed  CPU  Time. 

Reads  impact  parameters. 

Reads  laminate  materials  and  orientations.  Formu¬ 
lates  laminate  stiffness  matrices.  Computes  ef¬ 
fective  plane  stress  matrices  for  shell  element 
routines. 

Reads  modal  data.  Determines  equation  numbers. 
Sets  impacted  mode  numbers  to  degree  of  freedom 
numbers. 

Sets  storage  for  shell  element  routines. 

Supervises  shell  element  formulation. 

Shell  element  routines. 


Finds  global  stiffness  band  width  and  saves  elemental 
matrices  for  global  stiffness  formulation. 


ERROR 


INL 

ADDSTF 

ADSTF2 

STEP 

STATIC 

ADDMAS 

MASS IN 


SOLSTP 

TRIFAC 

REDVK 

PR 

STRREC 

PLYSTR 

FAIL 

BSTIF 

INVRTS 

MATPRD 

STOPPR 


Compares  required  storage  and  available  storage. 
Reads  and  processes  static  nodal  loads. 

Forms  global  stiffness,  mass  and  force  matrices. 
Computes  impact  velocity  for  momentum  conservation. 
Modifies  existing  global  stiffness  matrix  for 
damage  incorporation. 

Supervises  time  integration  analysis.  Sets  in¬ 
ternal  storage  for  integration  analysis. 

Supervises  solution  for  static  displacement  vector. 
Converts  blocked  mass  and  force  vectors  (from 
ADDSTF)  to  single,  unblocked  vectors. 

Reads  mass  vector  into  core  after  damaged  element 
reformulations.  Modifies  global  mass  and  stiffness 
for  impact  mass  separation. 

Performs  time  integration  analysis. 

Decomposes  stiffness  matrix. 

Solves  for  displacement  vector. 

Prints  displacements,  velocities  and  accelerations. 
Supervises  stress  recovery  procedure. 

Performs  composite  stress  analysis. 

Performs  failure  analysis. 

Computes  properties  for  interlaminar  shear  cal¬ 
culations  in  damaged  elements. 

Matrix  inversion  routine. 

Matrix  multiplication  routine. 

Terminates  program  execution. 


Ill 
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